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ABSTRACT 

Creep  crack  growth  in  a  nickel  base  alloy  at  elevated  temperatures 
was  analyzed  through  a  hybrid  experimental-numerical  (HEN)  procedure. 

This  HEN  procedure  consisted  of  simultaneous  use  of  creep  crack  growth 
test  displacement  data  from  center  cracked  plate  specimens  of  IN-100  at 
1350°  F  and  a  theoretical  finite  element  model  of  the  test  specimen. 

A  two-dimensional  (constant  strain  triangular  )  finite  element  program 
was  developed  which  accounts  for  both  nonlinear  viscoplastic  material 
behavior  and  changing  boundary  conditions  due  to  crack  growth.  Three 
viscoplastic  material  models  —  (1)  Malvern  Flow  Law,  (2)  Norton's  Creep 
Law,  and  (3)  Bodner-Partom  Flow  Law  —  were  incorporated  into  the  program. 
These  time  dependent  material  models  were  numerically  integrated  through 
time  by  a  linear  Euler  extrapolation  technique.  A  variable  time  step 
algorithm  was  included  that  maximized  time  step  size  during  the  analysis 
while  maintaining  good  accuracy.  This  program  was  used  as  the  plane  stress 
theoretical  model  for  the  HEN  procedure  to  analyze  sustained  load  creep 
crack  growth. 

A  method  for  getting  creep  crack  growth  behavior  solely  from  high 
resolution  displacement  measurements,  in  conjunction  with  a  cracked  specimen 
model  which  utilizes  realistic  constitutive  relationships,  has  been  developed. 
The  constitutive  law,  in  the  form  of  the  Bodner-Partom  material  model,  was 
especially  tailored  to  the  nickel-base  alloy  studied  which  displays  time 
dependent  nonlinear  inelastic  behavior  at  elevated  temperatures.  It  has 
been  demonstrated  that  the  technique  can  be  applied  where  crack  extension  is 
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very  small  and  could  not  otherwise  be  resolved  by  conventional  experimental 
crack  measuring  techniques.  This  method  provides  realistic  monotonically 
increasing  crack  growth  values.  The  predictions  agreed  to  within  10%  of 
post-test  measurements. 

Crack  growth  rate  and  crack  growth  criteria  were  studied.  Cracktip 
strain  and  crack  opening  displacement  were  studied  in  the  HEN  results  for 
a  unique  parameter  controlling  crack  growth.  Because  the  parameters  were 
not  independent  of  time  due  to  apparent  environmental  degradation,  it  became 
necessary  to  establish  an  empirical  criterion  for  crack  growth  based  on  the 
best  fit  of  HEN  results.  A  damage  accumulation  criterion  based  on  creep 
rupture  formulations  was  also  developed  and  applied  with  promising  results. 

Several  crack  growth  rate  criteria  were  investigated,  one  of  which 
is  the  stress  intensity  factor.  The  K  criterion  matched  fairly  well  with 
an  extrapolation  of  published  results.  But  two  other  criteria,  based  on 
the  C*  integral  and  load  point  displacement  rate  which  are  closely  related 
theoretically,  were  found  to  be  ineffective  as  crack  growth  rate  criteria. 
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I  INTRODUCTION 


BACKGROUND 

The  United  States  Air  Force  currently  places  strict  fracture  mechanics 
requirements  on  airframe  construction  and  maintenance  [1.2].  These 
requirements  involve  detection  of  flaws  by  periodic  nondestructive  in¬ 
spection  and  then  predicting  the  remaining  useful  life  of  the  part  through 
specified  fracture  mechanics  techniques.  Consequently,  if  an  airframe  part 
is  examined  and  found  to  have  no  flaws  that  can  grow  to  critical  size 
prior  to  the  next  periodic  examination,  it  can  be  returned  to  service. 

In  contrast  to  the  airframe,  low-cycle-fatigue  limited  jet  engine  parts 
may  be  retired  from  service  when  no  flaws  have  yet  been  found  in  them. 

This  situation  occurs  because  the  retirement  of  engine  disks  is  based  on  a 
"crack  initiation"  criterion.  Under  this  criterion  all  components  of  a 
given  population  are  considered  to  have  failed  as  soon  as  a  crack  of  some 
finite  size  e.g.,  .031  inches,  has  statistically  formed  in  the  member 
of  the  population  which  has  minimum  strength  properties  [3],  No  attempt 
is  made  to  utilize  the  additional  life  associated  with  the  remaining 
population  members  which  have  statistically  higher  properties  and  are 
therefore  expected  to  be  uncracked. 

From  a  safety  standpoint,  this  approach  has  been  generally  very  successful. 
However,  for  real  materials  and  real  design  situations,  lifetimes  based 
on  time  crack  initiation  of  the  minimum  member  tends  to  be  extremely  con¬ 
servative  for  a  component  population. 

It  has  been  estimated  that  replacement  costs  for  low-cycle-fatigue 
limited  jet  engine  disk  components  could  reach  the  $100,000,000  level  by  the 
1980  to  1985  time  period  [4],  A  significant  reduction  of  this  cost  could  be 
realized  if  a  procedure  was  developed  to  provide  accurate  assessment  of  useful 
residual  life  in  retired  engine  disks.  This  procedure  would  require  both 
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improved  inspection  and  fracture  mechanics  life  prediction  techniques. 

The  present  research  is  aimed  at  developing  more  accurate  fracture 
mechanics  life  prediction  capabilities.  Various  controlling  aspects  of 
fatigue  crack  growth  are  strain,  stress,  stress  intensity,  temperature,  load 
application  frequency  and  environment.  Speidel  [5]  provides  an  excellent 
discussion  of  the  relative  effects  of  each  of  the  above  aspects  on  fatigue 
crack  growth  rates  at  high  temperatures.  It  was  shown  that  at  each 
elevated  temperature  there  is  a  critical  frequency  below  which  the  crack 
growth  rate  is  creep  dependent  (i.e.  dependent  on  exposure  time  to  load)  and 
this  creep  dependency  increases  with  decreasing  frequency.  Also  the 
effects  of  an  aggressive  environment  is  another  time-dependent  phenomenon 
that  results  in  a  frequency  dependence  of  the  crack  growth  rate  which  is  very 
similar  to  that  brought  about  by  creep. 

Various  engine  missions  may  include  long  dwell  times  at  high  stress 
levels.  Crack  growth  controlling  aspects  may  then  be  reduced  to  stress  and 
strain  levels,  stress  intensity,  temperature  and  environment  since  load 
cycling  is  not  occurring  during  these  dwell  periods. 

Due  to  the  high  stress  concentration  in  the  vicinity  of  the  crack  tip, 
(i.e.,  infinite  stress  concentration  using  linear  elastic  fracture  mechanics) 
and  the  high  temperature  environment  for  an  engine  disk,  the  stress-strain 
relations  for  the  material  are  non-linear  and  time-dependent.  The  high  stress 
concentration  causes  the  material  to  yield  and  envelope  the  crack  tip  with  a 
plastic  zone.  Simultaneously,  the  high  temperature  allows  the  material 
under  load  to  flow  with  time,  the  phenomena  known  as  creep,  (i.e.,  increase 
strain  with  no  increase  in  stress).  Also  the  environment,  both  temperature 
and  atmosphere,  may  be  changing  the  ductility  of  the  material  thereby 
lowering  its  ability  to  strain  prior  to  fracture. 
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Creep  crack  growth  studies  of  metals  in  a  vacuum  indicate  that  removal 
of  air  from  the  crack  can  reduce  the  crack  growth  rates  by  two  orders  of 
magnitude  [5].  This  result  demonstrates  the  importance  of  recognizing 
environmental  effects  in  addition  to  what  may  be  called  mechanical  effects 
on  crack  growth.  Environmental  effects  may  be  thought  to  dictate  the 
critical  level  of  strain  at  fracture  whereas  mechanical  effects  determine 
the  rate  at  which  the  material  moves  to  the  critical  strain  levels. 

The  present  work  mainly  deals  with  the  mechanical  aspects  of  crack 
growth  under  fixed  load  in  materials  that  deform  with  time  (creep  crack  growth) . 
However,  since  the  experimental  data  used  herein,  came  from  elevated 
temperature  tests  in  laboratory  air,  some  environmental  effects  were  also 
considered . 

A  theoretical  model  for  creep  crack  growth  under  varying  and  fixed 
loads  must  be  able  to  account  for  the  changing  boundary  conditions 
associated  with  crack  growth.  In  the  event  of  total  unloading  and  reloading 
between  two  different  fixed  loads  the  possibility  of  crack  closure  and 
separation  needs  to  be  taken  into  account.  These  changing  boundary 
conditions  coupled  with  nonlinear  time  dependent  material  behavior  are  well 
suited  for  the  finite  element  method. 

APPROACH 

The  modeling  effort  considered  herein  involved  developing  a  two- 
dimensional  (plane  stress/plane  strain)  nonlinear,  time  dependent,  finite 
element  program  to  investigate  creep  crack  growth  under  constant  load.  The 
finite  element  analysis  incorporated  the  constant-strain  triangular  element. 

The  nonlinear  time-dependent  material  constitutive  model  took  the  form  of  the 
Bodner-Partom  viscoplastic  flow  law  [6, 7, 8,9].  This  flow  law  was  integrated 
through  time  by  an  Euler  extrapolation  scheme  [10]  and  incorporated  into  the 
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overall  finite  element  program  by  means  of  the  residual  force  method  [11]. 
Material  constants  for  the  Bodner  material  model  were  determined  by 
Stouffer  [12]  to  best  match  the  behavior  of  IN-100  [13],  a  current  jet 
engine  turbine  disk  alloy. 

Time  step  sizes  of  the  Euler  scheme  were  maximized  subject  to  specified 
amounts  of  change  in  stress  and  strain  over  a  given  time  step.  This  time  step 
maximization  scheme  provided  the  ability  to  transition  from  small  fraction  of 
a  second  time  steps  for  the  load  up  phase  to  large  time  steps  of  the  order  of 
minutes  or  even  hours  for  the  constant  load  creep  phase.  This  variable  time 
step  capability  was  a  necessity  to  make  a  numerical  study  of  creep  crack  growth 
computationally  feasible. 

Crack  growth  and  possible  crack  closure  during  unloading  was  accounted 
for  by  simple  modifications  to  the  structural  stiffness  matrix.  These  simple 
modifications  were  made  possible  by  choosing  an  iterative  Gauss-Seidel 
linear  equation  solver  [14]  which  requires  no  explicit  factorization  of  the 
stiffness  matrix.  Hence  between  time  steps,  pertinent  terms  of  the  stiffness 
matrix  could  be  easily  changed  to  account  for  crack  growth  and  the  general 
procedure  continued  without  costly  matrix  factorization  time  required. 

The  finite  element  program  which  includes  the  capability  of  accounting 
for  material  creep  behavior  and  crack  growth  was  used  to  study  the  creep 
crack  growth  in  a  center  cracked  plate  test  specimen.  Several  finite 
element  models  were  incorporated  to  study  different  initial  crack  lengths 
in  the  plate  geometry.  These  models  were  subjected  to  various  loads  that 
were  chosen  to  coincide  with  a  parallel  experimental  program  conducted  by 
W.  Sharpe  [15]. 

The  primary  objectives  in  the  present  research  were  to  determine  the 
actual  rate  of  creep  crack  growth  in  test  specimens  from  experimental  dis¬ 
placement  and  compliance  measurements  and  to  determine  the  most  reliable 
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criterion  for  predicting  creep  crack  growth  in  a  typical  jet  engine  turbine 
disk  alloy.  Specifically  the  interest  was  in  IN-100  at  1350°F.  These 
objectives  were  to  be  accomplished  by  the  so-called  hybrid  experimental - 
numerical  procedure  [16].  In  this  procedure,  crack  growth  rates  would  be 
estimated,  imposed  on  the  finite  element  simulation,  and  then  the  finite 
element  model  displacement  results  would  be  compared  with  experimental  data 
for  the  same  geometry  and  loading  conditions.  A  similar  method  is  to  allow 
the  crack  to  grow  sufficiently  so  that  predicted  crack  opening  displacement 
rates  from  the  finite  element  model  match  experimental  data  for  the  same 
geometry  and  loading  conditions.  After  good  correlation  of  displacements 
between  model  and  experiment  was  achieved,  fracture  criteria  were  sought 
out  from  the  calculated  parameters  in  the  finite  element  model  such  as 
stress,  strain,  and  displacement.  Criteria  were  sought  which  could  match 
the  now  determined  experimental  crack  growth  rates,  displacement,  or 
displacement  rates  over  a  range  of  crack  length  and  load  levels. 

Once  a  reliable  creep  crack  growth  criterion  is  determined  and  found 
to  be  independent  of  specimen  geometry  it  can  then  be  applied  to  an  actual 
turbine  disk  specimen.  With  the  determination  of  flaw  sizes  in  the  disk 
through  nondestructive  examination  techniques,  the  crack  growth  criterion 
could  then  be  used  to  predict  the  remaining  time  for  these  flaws  to  grow 
to  critical  dimensions  under  constant  load  applications. 
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II  LITERATURE  REVIEW 

Creep  crack  growth  in  general  may  be  thought  of  as  subcritical 
crack  growth  in  a  material  that  deforms  with  time  under  constant  external 
load.  This  time  dependent  deformation  or  creep  may  be  reversible 
(i.e.,  anelastic  creep)  or  it  may  be  permanent  (i.e.,  plastic  creep). 

At  elevated  temperature  metals  generally  exhibit  nonlinear  time  dependent 
deformation.  Under  uniaxial  tensile  loading,  the  strain  in  a  smooth  bar 
increases  with  time  until  failure  ultimately  occurs.  Based  on  similar 
response  of  many  materials,  researchers  have  subdivided  the  creep  curve 
into  three  regions  as  shown  in  Fig.  II-l.  After  the  initial  instantaneous 
strain  eq,  materials  often  undergo  a  period  of  transient  response  where 
the  strain  rate,  e,  decreases  with  time  to  a  minimum  steady-state  value 
that  persists  for  a  substantial  portion  of  the  materials  life.  These  two 
regions  are  referred  to  as  transient  or  primary  stage  and  steady-state 
or  secondary  stage  respectively.  Final  failure  with  rupture  life  tR  occurs 
soon  after  the  creep  rate  begins  to  increase  during  the  third  or  tertiary 
stage  of  creep.  A  common  empirical  relationship  between  creep  strain  rate 
and  stress  in  the  secondary  stage  of  creep  is  given  as: 

€  -  ( cr)  (ii-D 

where  a  is  the  uniaxial  stress,  e  is  the  creep  strain  rate,  and  y and  B 
are  empirical  constants  chosen  to  match  creep  test  results. 

Creep  crack  growth  has  been  studied  using  viscoelastic  [17-20], 
viscoelastic-plastic  [21]  and  plastic  creep  material  models  for  metals  as 
indicated  in  a  literature  review  by  Fu  [22].  The  viscoelastic  modeling, 
a  form  of  anelastic  behavior,  pertains  mainly  to  nonmetals  such  as 
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elastomers,  polymers,  and  solid  rocket  propellants.  The  material  of 
interest  in  this  investigation  is  a  current  jet  engine  turbine  disk 
alloy  known  as  IN-100  [13].  In  the  present  study  all  time  dependent 
deformation  of  IN-100  is  treated  as  permanent  and  consequently  any 
anelastic  behavior  is  considered  negligible.  However,  a  brief  review 
of  some  anelastic  (i.e.  viscoelastic)  creep  crack  growth  studies  will 
be  addressed  later. 

The  next  two  review  sections  will  consider  recent  creep  crack  growth 
studies  of  an  analytical  or  closed  form  nature  as  well  as  applications 
of  the  finite  element  method  to  creep  crack  growth. 

ANALYTICAL  CREEP  CRACK  GROWTH  STUDIES 

Knauss  [17]  analytically  modeled  steady  crack  growth  in  a  visco¬ 
elastic  sheet.  In  his  study  the  plastic  zone,  which  was  assumed  small, 
was  accounted  for  by  prescribed  fixed  and  finite  stress  distributions  in 
the  cracktip  region.  No  interaction  between  cracktip  and  the  far  field 
stresses  were  allowed.  This  means  the  stress  profile  ahead  of  the  moving 
cracktip  remained  constant  and  independent  of  the  far  field  stresses. 
Magnitudes  and  gradients  of  the  stress  in  this  cracktip  region  were  studied 
along  with  two  crack  growth  criteria.  The  two  criteria  were  maximum 
strain  and  a  maximum  strain  energy  criterion. 

Schapery  [18,  19,  20]  performed  a  viscoelastic  crack  growth  analysis 
similar  to  Knauss  but  placed  no  significant  restrictions  on  the  nature 
of  the  failing  material  at  the  cracktip.  It  could  be  highly  nonlinear 
and  rate  sensitive.  An  energy  criterion  for  failure  was  also  used  here. 

Wnuk  [21]  included  plasticity  with  viscoelasticity  for  his  quasi¬ 
static  extension  of  a  tensile  crack  analysis.  A  "final  stretch"  crack 
growth  criterion  was  proposed.  This  criterion  postulates  that  the  amount 
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of  deformation  which  occurs  within  the  cracktip  region  or  process  zone 


during  the  time  interval  just  prior  to  decohesion  of  this  zone  is  a 
material  constant.  In  contrast  to  the  maximum  strain  criterion,  the 
"final  stretch"  criterion  is  path-dependent. 

Fu's  review  [22]  of  quasi-static  crack  growth  in  metals  at  elevated 
temperature  presents  four  different  creep  crack  growth  rate  equations: 
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where  a  is  the  crack  growth  rate,  K  is  the  linear  elastic  stress  intensity 
factor,  is  the  net  section  stress,  (e.g.  load  divided  by  remaining 
uncracked  ligament  1b  a  center  cracked  plate  geometry),  y  in  general  is 
the  load-point  displacement  rate,  and  C*  is  a  line  integral  related  to 
the  rate  of  change  of  potential  energy  release  per  unit  of  crack  growth 
[22-25].  The  C*  integral  also  discussed  in  Appendix  C,  is  obtained 
directly  from  Rice's  J  integral  by  introducing  strain  rate  and  dis¬ 
placement  rate  instead  of  strain  and  displacement  such  that: 
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•  traction  vector 

u^,  ,  -  displacement  and  displacement  rate  respeclvely 

eij’  eij  "  8train  an<*  strain  rate  respectively 
I*,  x,  ds  -  see  Fig.  II-2 

There  are  available  experimental  data  supporting  any  of  the  rate 
Eqs  II-2  through  II-5  as  listed  by  Fu. 

The  following  conclusions  have  been  extracted  from  the  literature: 

1.  The  stress  exponent  B  in  Eq  II-l  plays  an  Important  role  in 

determining  the  characterization  of  the  crack  tip  behavior.  For 
8  5  the  stress  intensity  factor  approach,  Eq  II-2  may  be  used, 

and  for  8  _>  7  the  net  section  stress  approach  may  be  used  [26,27]. 

2.  Critical  test  conditions  for  evaluating  a  creep  crack  growth 
criterion  should  consist  of  at  least  two  geometries  which  have 
different  stress-intensity-factor  divided-by-stress  ratios  [23], 

Some  crack  growth  criteria  have  been  found  to  be  dependent  on 
geometry  and  therefore  have  no  general  application.  Use  of  two  or 
more  test  geometries  helps  determine  how  dependent  a  crack  growth 
criterion  is  on  geometry. 

3.  Creep  crack  growth  results  from  two  competing  processes.  These 
processes  are  (1)  growth  and  coalescence  of  defects  which  contribute 
to  crack  advancement  and  (2)  the  creep  deformation  process  that  causes 
retardation  and  even  arrest  of  crack  growth  [25]. 

A.  Creep  crack  growth  rates  are  very  sensitive  to  environmental 
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effects.  Removal  of  air  from  the  crack  can  reduce  the  crack  growth 
rates  by  two  orders  of  magnitude  [5]. 

5.  Crack  opening  displacement  crack  growth  theories  indicate 
that  failure  times  due  to  creep  crack  growth  are  controlled  by 
the  stress  intensity  factor  at  large  stresses  and  by  net  section 
stresses  at  very  low  stresses  [28].  However,  a  counter  viewpoint 
is  stated  in  [29]  where  it  is  concluded  that  creep  crack  growth 

does  not  correlate  well  with  the  stress  intensity  factor  at  relatively 
high  stress  levels. 

6.  Creep  crack  growth  rates  correlate  with  the  energy  rate  integral  C*. 
This  method  holds  great  'promise  for  design  calculations  because 

C*  can  be  calculated  using  finite  element  analysis  as  well  as 
measured  empirically  in  constant  displacement  rate  tests  [24]. 

7.  An  approximation  to  C*  in  the  form  of  the  product  of  net 
section  stress  and  load  line  displacement  rate,  referred  to  as  J, 
gives  good  correlation  of  creep  crack  growth  rates  in  specimens 
of  different  geometries  [29]. 

8.  Crack  growth  theories  generally  fall  into  one  of  two  categories. 
Either  they  are  of  an  energy  nature  (e.g.  J  or  C*  integrals), 

or  they  deal  with  some  localized  crack  tip  parameter  such  as  strain 
or  crack  opening  displacement. 

FINITE  ELEMENT  ANALYSIS  RELATED  TO  QUASI-STATIC  TIME  DEPENDENT  CRACK  GROWTH 
The  general  technique  of  approximating  a  continuum  with  simple 
discrete  elements  such  as  uniaxial  bar  elements  dates  back  to  the  1940's 
[30].  The  history  of  the  finite  element  method  in  structural  analysis, 
as  it  is  known  today,  is  well  described  by  Zienkiewicz  [31].  The 
application  of  the  method  to  the  problem  of  nonlinear  material  behavior 
has  also  been  developed  [31].  The  purpose  of  this  section  is  to  briefly 


review  the  use  of  the  finite  element  method  for  the  stress  analysis  of 
cracked  plates  where  nonlinear  time  independent  and  time  dependent 
materials  models  were  employed. 

ELASTIC- PLASTIC  (TIME  INDEPENDENT)  ANALYSIS.  The  finite  element 
method  has  been  widely  used  to  determine  the  stress  and  strain  fields 
around  cracks  in  nonlinear  materials  where  time  Independent  elastic- 
plastic  materials  models  were  Incorporated  [32-45].  Several  of  these  use 
the  "Initial  stress"  or  "initial  strain"  approach  to  elastic-plastic 
modeling  as  described  in  Appendix  A.  During  the  last  few  yearB  many 
reports  have  been  published  on  this  subject.  This  review  will  not  cover 
all  these  publications,  however,  it  will  concentrate  on  two  areas;  (1) 
crack  tip  element  selection  and  (2)  crack  growth  modeling  by  the  finite 
element  method. 

There  are  several  choices  to  consider  when  selecting  finite 
elements  to  model  the  vicinity  of  a  crack  tip.  These  choices  might  best 
be  classified  into  two  categories.  The  first  category  is  to  use  the 
same  element  type  incorporated  in  the  model  remote  from  the  crack  tip 
(e.g.  constant  strain  triangles).  The  second  category  is  to  use  a 
special  crack  tip  element  that  has  functions  to  include  the  specific 
stress  or  strain  singularity  desired  at  the  crack  tip. 

There  are  of  course  benefits  and  disadvantages  to  each  choice. 

Using  the  same  element  such  as  a  constant  strain  triangle  both  at  the 
crack  tip  and  remote  definitely  has  the  benefit  of  simplicity.  However, 
the  cost  of  this  simplicity  is  the  requirement  to  use  large  numbers  of 
such  elements  around  the  crack  tip  to  obtain  acceptable  results.  Also 
in  a  fully  plastic  material,  some  elements  do  not  accurately  model 
incompressible  strain  behavior  [46].  The  bilinear  displacement 
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quadralateral  Is  most  susceptible  to  this  inaccuracy,  furthermore  the 
constant  strain  triangle  is  one  of  the  least  susceptible. 

The  selection  of  a  special  cracktlp  element  requires  knowledge 
of  the  stress  or  strain  singularity  around  the  cracktlp.  Several  authors 
have  developed  the  inverse  square  root  "r"  singularity  for  linear 
elastic  cracktlp  behavior.  Accordingly  many  special  elements  have  been 
addressed  to  this  singularity.  In  a  paper  entitled  "Cracktlp  Finite 
Elements  Are  Unnecessary"  [A7],  the  authors  describe  the  modification  of 
the  eight  noded  isoparametric  element  such  that  it  Incorporated  the 
1/fr  stress  singularity.  Therefore  any  existing  finite  element  program 
that  had  the  eight  noded  isoparametric  element  in  it  also  effectively  has 
one  form  of  a  special  cracktlp  element  capability  for  elastic  analysis. 

Further  considerations  of  the  local  cracktlp  stress  and  strain 
distributions  in  an  elastic-plastic  material  have  led  to  the  expressions 
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where  K^,  K£  are  plastic-intensity  factors  and  n  is  a  strain  hardening 
exponent.  This  type  of  singularity  was  embedded  in  a  cracktlp  element 
and  used  to  evaluate  the  plastic  intensity  factors  as  a  function  of 
applied  loading  [49].  The  advantage  of  this  method  was  that  the  large 
strain  gradients  in  the  cracktlp  region  are  accounted  for  by  Eq  11-10. 
For  elastic-plastic  and  creep  type  materials  the  singularity  changes 
between  initiation  and  growth  of  the  crack  and  in  general  is  not  known. 
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The  development  of  a  special  crack  tip  element  in  the  case  where  the 
singularity  is  not  known  to  begin  with  would  be  an  extensive  undertaking 
by  itself  [37]  without  bringing  in  the  additional  complexity  of  a 
moving  crack  tip. 

Finite  element  researchers  have  considered  the  crack  growth  problem. 
Kobayashi,  Chiu,  and  Beeukes  [43]  analyzed  an  extending  crack  under 
monotonically  increasing  load.  Crack  extension  was  achieved  by  applying 
a  relief  force  equal  in  magnitude  but  opposite  in  direction  to  the 
restraining  force  at  the  crack  tip  node.  This  relief  force  was  applied 
in  100  equal  increments  or  in  one  single  increment.  The  crack  opening 
displacements  at  the  node  adjacent  to  the  crack  tip  computed  by  the 
single  increment  method  were  less  than  5  per  cent  smaller  than  the 
corresponding  displacements  computed  with  the  100  increment  method. 

Lee  and  Liebowitz  [35]  ising  a  similar  technique  demonstrated  that  plastic 
strain  energy  increased  linearly  with  crack  length  as  the  crack  grew 
in  their  analysis.  Anderson  [36]  also  made  use  of  relaxing  the  crack  tip 
node  force  incrementally  to  simulate  crack  growth  in  the  finite  element  model. 

Shih,  DeLorenzi,  and  Andrews  [42]  analyzed  crack  initiation  and 
stable  crack  growth  in  elastic-plastic  material  by  special  use  of  eight- 
node  isoparametric  elements.  Initial  crack  tip  blunting  was  modeled 
followed  by  crack  extension.  Crack  extension  was  modeled  by  a  combination 
of  crack  tip  node  shifting  and  then  releasing  it  to  move  on  to  the  next 
element's  nodes.  Moving  these  nodes  within  elements  at  the  crack  tip 
required  approximations  to  be  made  about  Gauss  point  location  and  relocation. 
These  approximations  are  not  necessary  when  only  a  node  release  method 
and  elements  such  as  constant  strain  triangles  are  employed  to  model 
crack  growth.  The  1 Jr  strain  singularity  provided  by  the  special  use 
of  these  eight  noded  isoparametric  elements  may  better  model  the  crack  tip 
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singularity  for  a  fixed  crack  length  and  a  theoretical  continuum. 

However,  when  the  crack  initiates  and  grows  in  a  creep  type  material 
the  strain  singularity  is  unknown,  especially  when  considering  a  grain 
structure  around  the  crack  tip  rather  than  a  continuum  and  the  fact  that 
a  creep  crack  follows  an  intergranular  path. 

VISCO-PLASTIC  AND  CREEP  (TIME  DEPENDENT)  ANALYSIS.  Several  references 
have  been  found  on  viscous  or  time  dependent  material  models  being 
incorporated  into  finite  element  program  [50-60],  In  general  these 
material  models  may  be  similar  to  Eq.  II-l  for  pure  creep  with  the 
addition  of  time  independent  elastic-plastic  relationships,  or  they 
may  have  short  term  response  viscoplastic  relationships  that  only  model 
the  load  up  phase.  Zienkiewicz  [10]  proposed  placing  a  short  term 
response  viscoplastic  equation  in  series  with  a  long  term  response  creep 
law,  similar  to  Ea.  II-l.  This  would  be  a  unified  time  dependent  material 
model  where  no  direct  coupling  is  assumed  between  short  term  plastic 
strains  and  long  term  creep  strains. 

Only  a  few  papers  have  been  found  to  date  that  use  a  time  dependent 
material  model  and  the  finite  element  technique  to  model  crack  growth. 

Ohtani  and  NaKamara  [61]  analyzed  creep  crack  propagation  with  an  elastic¬ 
secondary  creep  material  model.  The  secondary  creep  law  was  identical 
to  Eq .  II-l  in  its  uniaxial  form.  A  critical  strain  criterion  (i.e.  average 
effective  plastic  strain  at  the  crack  tip)  was  used  for  determining  the 
time  to  grow  the  crack  in  the  finite  element  model  of  a  center  cracked 
plate.  Crack  tip  nodes  were  released  to  simulate  crack  growth.  Goodall 
and  Chubb  did  a  similar  analysis  on  the  compact  tension  specimen  [62). 

A  critical  strain  crack  growth  criterion  was  again  used.  Finally 
Zaphir  and  Bodner  [63]  incorporated  the  time  dependent  Bodner-Partom 
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viscoplastic  flow  law  into  the  NONSAP  finite  element  code  to  analyze 
a  double-centilever-beam  cracked  geometry.  In  this  case,  high  loading 
rates  were  studied  over  short  time  periods  relative  to  creep  analyses. 
Consequently  the  "recovery  term",  as  described  in  a  later  section  and 
used  to  best  model  creep,  was  not  Included.  No  crack  growth  was  allowed 
here. 

In  each  of  the  finite  element  solutions  referred  to  previously, 
a  unified  time  dependent  material  model  that,  not  only  accurately 
models  the  short  term  load  up  stage  of  material  response,  but  also 
smoothly  transitions  into  the  pure  creep  stage  was  never  considered. 

In  addition,  the  hybrid  experimental-numerical  technique  was  not 
used  with  high  resolution  experimental  crack  opening  displacement  data. 

It  was  anticipated,  in  the  present  research,  that  the  combination  of  a  more 
realistic  time  dependent  material  model  and  high  resolution  test  data 
would  result  in  a  much  better  understanding  of  what  controls  creep  crack 
growth  than  provided  by  these  prior  analyses. 
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Ill  FINITE  ELEMENT  COMPUTER  PROGRAM  DEVELOPMENT 


This  section  describes  the  development  and  validation  of  a  finite 
element  computer  program  called  "VISCO".  VISCO  is  a  two  dimensional 
plane  stress/plane  strain  program  that  incorporates  three  different 
nonlinear  time  dependent  viscoplastic  material  models.  It  uses  constant 
strain  triangular  elements  with  the  option  to  release  fixed  nodes  and  thus 
has  the  capability  to  simulate  crack  growth.  Results  from  VISCO  are 
compared  with  other  published  solutions  to  check  its  validity. 

The  approach  selected  for  elastic-viscoplastic  analysis  with  the 
finite  element  technique  employs  the  "residual  force  method"  [11], 

(see  Appendix  A  for  a  complete  discussion  on  this  particular  approach) . 

In  the  residual  force  method  the  elastic  stiffness  is  used  during  the  entire 
analysis  and  any  nonlinear  elastic-viscoplastic  deformation  that  occurs 
must  be  accounted  for  by  developing  so-called  plastic-load  vectors  that 
add  to  the  force  side  of  the  governing  equilibrium  equation.  In  general, 
the  matrix  equation  which  governs  the  response  of  a  discretized  structure 
can  be  written  as 

lk]  *pr  +  ( ii i— i, 

where  [K]  is  the  elastic  stiffness  matrix  {l1  J1  is  the  generalized  nodal 
displacement  vector  for  the  ith  time  step,  (P  l1  is  the  load  vector  after 
the  ith  time  step  due  to  external  forces,  and  {Q  ^  is  the  plastic-load 
vector  computed  from  plastic  strains  accumulated  prior  to  the  ith  time  step. 
For  each  element,  these  plastic-load  vectors  are  self-equilibrating.  The 
viscoplastic  strain  rate  expressions  which  develop  plastic  strains  with 
time  under  load  are  described  below. 
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VISCOPLASTIC  MATERIALS  MODELS 


In  solid  mechanics  it  is  customary  to  separate  the  two  important 
groups  of  phenomena  described  respectively  by  "creep"  and  "plasticity". 

The  first  includes  all  time  dependent  effects  and  results  in  creep  strains 
accumulating  at  a  finite  rate.  The  second  group  develops  permanent 
(plastic)  strains  instantaneously,  since  time  does  not  enter  directly 
into  consideration  as  in  the  elastic-plastic  approaches  given  in  Appendix  A. 
From  a  physical  point  of  view,  creep  and  plasticity  cannot  be  treated 
separately  as  only  the  combined  effect  is  measurable.  Also  the  concept 
of  time  independent  or  instantaneous  plasticity  is  only  a  convenient 
mathematical  approximation  and  not  experimentally  based. 

Viscoplasticity,  in  a  complete  sense,  is  the  combination  of  these  two 
strain  groups  into  a  unified  plastic  strain  rate  model.  A  model  with  this 
capability  is  the  Bodner-Partom  viscoplastic  flow  law  [6-9]  from  hereon 
referred  to  as  the  Bodner  model.  The  superposition  of  Malvern's  overstress 
law  [66,67]  with  Norton's  law  for  secondary  creep  has  also  been  proposed 
as  a  unified  viscoplastic  flow  law  [10], 

Each  of  these  flow  laws  has  been  incorporated  into  VISCO  by  assuming 
small  strains  and  decomposing  the  total  strain  rate  into  elastic  (reversible) 
and  plastic  (nonreversible)  components. 

£.-j  r  +  cni-2) 

which  in  general  are  both  non-zero  for  all  loading/unloading  conditions. 
Anelastic  stresses  and  strains  corresponding  to  time  dependent  reversible 
deformations  with  energy  losses  are  not  considered  in  this  formulation 
and  are  assumed  to  be  relatively  unimportant.  This  assumption  seems  quite 
justified  based  on  the  good  correlation  between  predicted  and  test  results 
to  be  shown  later. 


IP 


The  elastic  strain  rate  £  is  related  to  the  stress  rate  by  the  time 
derivative  of  Hooke’s  law.  The  plastic  strain  6  assuming  incompressibility 
and  isotropy,  is  taken  to  follow  the  Prandtl-Reuss  flow  law  of  classical 
plasticity 

£.?  =  2  ( III— 3) 

where  S^.  are  the  components  of  the  deviatoric  stress  tensor  and  X  is  a 
scalar  that  reflects  the  viscosity  of  the  material.  The  specific  form  of  X 
is  presented  below  for  each  of  the  viscoplastic  flow  laws. 

MALVERN  (OVERSTRESS)  FLOW  LAW.  A  portion  of  a  total  viscoplastic 
model  that  accounts  for  the  so-called  instantaneous  plasticity  during 
loading  might  take  the  form  as  given  by  Malvern  [66,67],  otherwise  known 
as  the  "overstress"  model: 


.  P 


0  I*  cr«  -  &(€e) 


(III-4 ) 


where  is  a  fluidity  constant  whose  magnitude  will  determine  the  strain 
rate  sensitivity  of  the  model,  see  Fig  III-l,  is  the  effective  stress 
defined  as  3  where  Jj  the  second  invariant  of  the  deviatoric  stress 
defined  as  =  l/2S^S^j  ,  and  is  the  strain  hardening  yield  stress 

shown  to  be  a  function  of  the  effective  plastic  strain,  e^,  defined  incrementally 
as  de^  d£^~def^ .  The  strain  hardening  stress  function,  within  VISCO, 
takes  one  of  two  forms,  either  a  linear  relationship  such  as 

<r(€e)=  ^  +  hi'*  (IH-5) 


or  a  Ramberg-Osgood  type 
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Figure  III-l.  Variation  of  Malvern  Model  Behavior  With 

"Yp  Under  Constant  Strain  Rate  Loading  Conditions 
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if  K(i ')-£  5-. 
if  K(4)~  <  %  nn-6) 

where  ao  is  the  initial  yield  stress,  the  constant  H*  represents  the  slope 

of  the  stress  versus  plastic  strain  curve,  and  R  and  m  are  constants  of 

a  Ramberg-Osgood  type  stress-strain  curve.  In  EqIII-6,  note  that  if  the 

effective  plastic  strain  is  at  or  near  zero  such  that  the  function 

e 

R(e^)  mis  less  than  the  initial  yield  stress,  crQ,  then  o(e^)  is  set 

equal  to  a  . 

o 

Implementation  of  the  Malvern  model  then  requires  selection  of 
Eq  III-5  or  Eq  III-6.  This  selection  would  depend  on  the  best  fit 
of  the  material's  uniaxial  stress-strain  curve  developed  under  strain  rates 
at  or  near  the  lowest  strain  rate  expected  to  be  modeled  with  the  Malvern  Law. 
If  Eq  III-5  were  chosen,  the  initial  yield  stress,  oq,  and  slope,  H’  would 
be  determined  from  this  experimental  curve.  A  similar  determination  would 
be  done  if  Eq  III-6  were  selected.  The  fluidity  constant,  Yp ,  would 
be  chosen  to  best  reflect  the  strain  rate  sensitivity  of  the  material 
(see  Fig.  III-l)  displayed  by  experimental  stress-strain  data  developed 
at  high  strain  rates. 

The  Malvern  model  may  also  be  used  to  perform  time  independent  elastic- 
plastic  solutions.  In  this  case  time  becomes  a  fictitious  parameter  and 
thus  the  fluidity  constant,  Yp,  may  take  on  any  non-zero  positive  value. 

The  elastic-plastic  solution  is  the  steady  state  value  of  the  stresses, 
strains,  and  displacements  after  the  load  has  been  applied.  This  has  been 
found  to  occur  in  approximately  30  timesteps  after  maximum  load  application 
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unless  total  section  yielding  develops. 

An  Euler  linear  extrapolation  scheme  is  employed  within  VISCO  for  the 
time  integration  of  the  viscoplastic  strain  rate  expressions.  The 
Malvern  model  is  integrated  in  VISCO  as  follows 


i  f  Ce"1  <  O'  ) ' "  * 
if  o^'1  -r 


( I II— 7 ) 


(1II-8) 


(111-9) 


o;  +  H'*  C^c  ) 

°r  .  X 

Rt<0';r 


(III-10) 


where  the  superscript  i  refers  to  the  timestep  and  a  subscript  i  refers 
to  specific  components  of  stress  or  strain. 

NORTON'S  LAN  FOR  SECONDARY  CREEP.  Another  portion  of  a  total  visco¬ 
plastic  flow  law  that  complements  the  Malvern  model  and  accounts  for  long 
term  creep  is  given  by  Norton's  creep  law  [68]  written  in  multiaxial  form  as 


Oe)  T 


(iii-id 


where  y and  8  are  constants  determined  from  uniaxial  creep  test  results. 
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Creep  test  data  at  two  different  stress  levels  are  required.  A  straight 
line  is  fitted  to  each  test's  secondary  stage  of  creep  strain  plotted  versus 
time,  see  Fig  II-l.  The  slope  of  this  line  or  strain  rate  and  stress 
level  from  each  test  is  substituted  into  Eq  III-ll  which  provides  two 
equations  to  solve  for  the  two  unknown  constants.  After  taking  the 
natural  logarithm  of  eq  III-ll,  the  simultaneous  solution  of  the  two 
equations  is  straight  forward. 

The  Euler  extrapolation  scheme  for  integrating  Eq  III-ll  in  VISCO 
is  similar  to  the  one  given  for  the  Malvern  model  and  will  not  be  re¬ 
peated  here. 

It  would  appear  that  the  combination  of  the  Malvern  and  Norton  laws 
would  provide  a  complete  viscoplastic  flow  law.  However,  metals  at 
elevated  temperature  have  been  observed  to  display  a  phenomenon  called 
"recovery"  [69].  Recovery  is  the  softening  of  cold-worked  metal  or  it 
may  be  characterized  as  a  fading  memory  of  prior  strain  hardening. 

In  creep  crack  growth  at  elevated  temperature  and  under  constant 
load,  consider  material  well  ahead  of,  but  in  the  path  of,  the  crack. 

During  load  up  this  material  will  plastically  deform  and  strain  harden 
depending  on  its  proximity  to  the  initial  crack  tip.  Subsequently,  during 
the  sustained  load  phase,  the  phenomena  of  recovery  will  allow  this  material 
to  soften  prior  to  the  arrival  of  the  crack  tip  at  which  time  additional 
loading  will  occur.  The  amount  of  recovery  prior  to  the  arrival  of  the 
crack  tip  should  then  have  some  effect  on  the  values  of  stress  and  strain 
developed  around  the  crack  tip  when  it  reaches  the  material  being  considered 
Therefore,  this  investigation  will  focus  on  the  following  viscoplastic  flow 
law  which  does  include  the  phenomenon  of  recovery. 

BODNER  VISCOPLASTIC  FLOW  LAW.  In  this  formulation  by  Bodner  and 


second  invariants  by  making  use  of  the  square  of  Eq.  III-3. 


(III-12) 


where  is  the  second  invariant  of  the  plastic  strain  rate  and  is  the 
second  invariant  of  the  deviatoric  stress. 

Rather  than  specifying  a  specific  yield  criterion  as  in  classical 
plasticity,  this  formulation  by  Bodner  and  Partom  is  based  on  the  assumption 
that  a  continuous  functional  relationship  exists  between  the  plastic 
deformation  rate  and  the  stress  invariants,  i.e. 


■>  hk  i T) 


(III-13) 


where  Z^  are  one  or  more  internal  (viscoplastic)  state  variables  and  T  is 
the  temperature.  Introducing  Eq.  III-13  into  Eq .  III-12  and  solving 
for  A  gives 


A  =  [Dr(^,*KnT>/T2  V 


(III-14) 


The  general  form  for  the  evolution  equation,  i.e.  history  dependence, 
of  the  viscoplastic  state  variables  Z^  is 


Z  =  Fk 


(III-13) 


For  conditions  of  uniaxial  stress  of  constant  sign,  the  hardened 
state  with  respect  to  plastic  flow  is  assumed  to  be  represented  by  a 
single  state  variable  Z  which  depends  on  plastic  work.  This  single  state 
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variable  Z  also  corresponds  to  isotropic  hardening.  Additional  state 


variables  are  necessary  for  such  characteristics  as  kinematic  hardening 
which  will  not  be  employed  here. 

The  particular  form  chosen  for  (^.Z.T)  was  motivated  by  the 
equations  of  dislocation  dynamics  and  given  by  Bodner  and  Partom  as 


(III-16) 


The  factor  (n  +  l)/n  was  introduced  at  an  early  stage  in  the  development 

of  the  equations  for  numerical  purposes  and  only  affects  the  numerical 

values  of  some  of  the  material  constants.  The  constant  D  is  described 

o 

as  the  limiting  value  of  the  plastic  strain  rate  in  shear.  Its  value  can 

be  arbitrarily  chosen  and  is  usually  taken  to  be  the  same  large  number 

4  _i 

for  all  materials,  A  value  Dq  =  10  sec  is  generally  adequate  except 
for  conditions  of  very  high  rates  of  straining. 

The  parameter  n  controls  strain  rate  sensitivity  and  also  influences 
the  overall  inelastic  level  of  the  stress-strain  curves.  A  decrease 
of  n  leads  to  increasing  strain  rate  sensitivity  and  lowering  of  the 
level  of  the  stress-strain  curves. 

As  stated  previously,  Z  is  assumed  to  be  a  function  of  plastic  work, 
Wp,  and  the  following  relationship  is  introduced 

Z  r  l( Wp)  r  2,  -  Ze)  e*p[-m  Wp]  (III-17 


The 


quantity  Z^  is  the  maximum  value  of  Z  which  is  necessary  if  the 


deformations  do  not  revert  to  elastic  behavior  at  large  values  of  W  . 

P 


Z  is  the  value  of  Z  for  which  W 
o  p 


o  and  can  therefore  be  the  initial 
state  point  from  which  plastic  work  is  measured.  It  is  noted  that  the 


general  function,  Eq  III-17,  would  be  a  basic  material  property  and 


that  Wp  is  the  relative  amount  of  plastic  work  done  from  some  initial 
state,  (i.e.  W  is  not  an  absolute  parameter). 


Wp  -  J  S,j  <|t 


(III-18) 


The  quantity  m  in  Eq  III-17  is  a  material  constant  that  relates  to 
the  rate  of  work  hardening. 

At  high  temperatures,  it  is  generally  necessary  to  consider  the 

thermal  recovery  of  hardening  generated  by  plastic  deformation.  In 

this  case  the  plastic  work,  W  ,  is  redefined  as  follows: 
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=  fe,  €P 
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where 
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(III-19) 


(III-20) 


where  is  the  stable,  (i.e.  non-work  hardened)  value  of  Z  at  a  given 

temperature,  and  A  and  r  are  additional  material  constants  chosen  to 

match  the  models  behavior  to  creep  test  data  as  was  done  for  the  Norton 

model.  Note  that  the  second  term  on  the  right  of  Eq  III-19  (i.e.  the 

recovery  term)  makes  a  negative  contribution  to  W  due  to  the  negative 

P 

sign  on  A,  since  Z  is  always  greater  than  or  equal  to  Z^ . 

The  recovery  term  in  Eq  III-19  is  essential  if  the  material  model  is 
to  be  able  to  represent  secondary  creep.  Secondary  creep  is  the  balanced 


condition  when  the  rate  of  work  hardening  equals  the  rate  of  thermal 
recovery  or,  setting  the  time  derivative  of  Eq  III-19  to  zero, 

VVp  =  S.  j  €?■  +  r  o  (III-21) 

At  relatively  high  strain  rates,  the  thermal  recovery  in  Eq  III-19 
is  relatively  unimportant  and  the  steady  state  condition  is  realized 
when  Z  reaches  its  saturation  value  Z^. 

Again,  VISCO  employs  the  Euler  extrapolation  scheme  for  the  numerical 
time  integration  of  the  Bodner  equations.  During  each  time  step  Eq  III-3 
and  Eqs  III-14  through  III-20  are  integrated  as  follows  for  each  element 


Z '  r  2,-  (vOe,lP[‘m  Wp  J 

(III-22) 

(°*)=  '*?[-  (  3  Z--<) 

(III-23) 

'  lo:>!/ i r' 

(III-24) 

fJ  «.!}■'= 

(III-25) 

i;.c  -  -  A(^-^-)rZ, 

(III-26) 

Wf'  =  V/p^+JS,^’  *0  J 

(III-27) 

where  the  superscript  i  refers  to  the  timestep  and  a  subscript  i  refers 
to  a  specific  material  constant,  Z^,  or  to  specific  components  of  stress 
and  strain. 


The  specific  material  constants  required  for  the  Bodner  model  in  this 
investigation  were  determined  to  best  fit  the  behavior  of  IN-100  at  the 
temperature  of  1350°F.  This  material  was  characterized  by  performing 
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uniaxial  tensile  stress-strain  tests  and  creep  tests  at  different  stress 
levels  at  1350°F,  The  following  constants  for  the  Bodner  model  were 
developed  by  Stouffer  [12]  and  the  details  of  this  procedure  are  summarized 
in  Appendix  0.  The  A  and  r  constants  defined  by  Stouffer  are  different 

than  the  constants  used  in  the  present  formulation: 

4  -1 

D  =  10  sec 
o 

n  =  3.50 

Z  =  224.4  ksi 
o 

Z  =  251.5  ksi 
m  =  3. 750  ksi  ^ 

This  first  group  of  constants  is  based  on  stress-strain  curve  data. 

A  =  1.142  x  10  2  sec  ^ 
r  =  3.52 
Z±  =  100  ksi 

This  second  group  of  constants  is  based  on  creep  test  data.  The  elastic 
modulus  at  1350°F  was  determined  to  be 

E  «  26.3  x  10  ^  psi 

and  Poisson's  ratio  was  arbitrarily  chosen  as 

v  =  .  3 

Figure  III-2  shows  a  comparison  of  the  stress-strain  behavior  of 

the  Bodner  model  (using  the  given  IN-100  material  constants)  and  experimental 

stress-strain  data.  Each  of  the  curves  displays  the  response  of  the 

Bodner  model  under  the  loading  condition  of  a  given  constant  strain  rate. 

-3  -1 

The  experimental  data  were  generated  at  a  strain  rate  of  .83  x  10  sec 
and  compare  quite  well  with  the  Bodner  curve  for  the  same  strain  rate. 

Some  variation  occurs  in  the  region  of  initial  inelastic  behavior  but 

this  difference  Is  small. 
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Figure  III-2.  Stress-Strain  Behavior  of  Bodner  Model 
(IN-100  Constants) 


Figures  1 1 1  —  3  and  1 1 1—4  compare  creep  behavior  of  the  Bodner  model 
with  experimental  data  at  the  stress  levels  of  127.3  ksi  and  72  ksi, 
respectively.  Note  the  initial  experimental  curve's  slope  or  strain  rate 
is  duplicated  by  the  Bodner  model  in  both  figures.  However,  the  strain 
magnitudes  differ  somewhat  due  fn  the  initial  time  required  for  the 
Bodner  model  to  reach  steady  state  creep  at  these  stress  levels.  Also 
as  the  experimental  strain  rate  increases  with  time  the  Bodner  model 
cannot  closely  follow  since  its  formulation  restricts  it  to  only  secondary 
type  creep  in  this  situatijn  (i.e.  constant  strain  rate  for  constant  stress), 
SOLUTION  PROCEDURE  FOR  ELASTIC-VISCOPLASTIC  STRUCTURES 

In  elastic-plastic  analysis,  it  is  necessary  to  apply  loads  incrementally 
to  satisfy  the  appropriate  yield  condition  (e.g.  von  Mises)  and  flow 
rule  (e.g.  Prandtl-Reuss)  associated  with  incremental  plasticity  [64], 
Similarly,  with  elastic-viscoplastic  behavior  an  incremental  procedure 
•is  required,  but  here  time  is  incremented  directly  while  load,  strain, 
stress,  etc.  are  incremented  indirectly  through  the  time  integration 


procedure.  The  algorithm  used  for  a  typical  time  step  in  the  elastic- 
viscoplastic  residual  force  method  [10]  is  summarized  as  follows: 

1.  Add  time  increment  dt*  to  the  preceding  time  t*  ^  to  obtain 
the  current  time  t  . 


2.  Compute  increments  of  plastic  strain,  {de^}1  =  {e^l^dt*  and 


add  to  preceding  plastic  strain,  {e?  } 1  = 


(c;,)1-1  +  IdcU)1.  In 
IJ  ij 


general  the  plastic  strain  rate,  {£]].},  is  a  function  of  the  material's 


viscosity  and  the  given  stress  level  (see  the  Viscoplastic  Material 


Models  section) . 

3.  Compute  the  plastic  load  vector  {Q}i  ^  =  /  ,[B][D]  (e^.i^dvol. 

vol  ij 


% 
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STRAIN 


Figure  III-3.  Creep  Behavior  of  Bodner  Model  At  127.3  Ksi 
(IN-100  Constants) 
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IN-100  Creeo(a)72.  ksi 
1 350° F,  Ref.  12 


345  min. 


4.  Compute  the  current  external  load  vector  {P } 1  =  {Pj^dt1  +  {P}*  ^ . 

5.  Compute  the  nodal  disp’  -cements  { U } 1  from  Eq.  III-l, 

{u}1  =  [k]“1£{p}1  +  {Q}1-1/. 

6.  Compute  the  current  total  strain  from  the  strain  dis¬ 

placement  relationship,  {e^}1  =  [BjfU}1. 

7.  Compute  the  current  stress  {o}^  as  follows,  j }  ^  -  {e^}^. 

8.  Check  the  timestep  size  in  terms  of  prescribed  stress  and  strain 
change  tolerances  per  time  step  (see  the  following  section  on 
Variable  Time  Step  Integration  of  Viscoplastic  Flow  Laws).  If 
these  tolerances  are  not  exceeded  the  time  step  size  may  be  increased 
for  the  next  time  step  or  left  the  same  value.  However,  if  the 
tolerances  are  exceeded,  the  time  step  size  is  reduced  and  steps 

1  through  8  are  repeated  for  this  same  time  step  in  an  effort  to 
satisfy  the  stress  and  strain  change  tolerances. 

9.  Repeat  steps  1  through  8  until  the  desired  simulation  time 
is  reached. 


VARIABLE  TIKE  STEP  INTEGRATION  OF  VISCOPLASTIC  FLOW  LAWS 

One  of  the  final  objectives  in  the  development  of  VISCO  was  to  be 
able  to  accurately  model  nonlinear  material  behavior  both  in  the  high 
strain  rate  load  up  stage  and  during  the  low  strain  rate  constant  load 
creep  stage.  For  stable  accurate  results,  the  time  step  size  must  be 
orders  of  magnitude  less  during  the  load  up  stage  compared  to  the  creep 
stage.  Therefore  to  be  computationally  feasible  in  a  large  problem  (many  degrees 
of  freedom)  some  method  is  necessary  to  determine  the  maximum  time 

step  during  each  stage  of  the  analysis  while  maintaining  reasonable 
accuracy.  Cormeau  [65]  investigated  the  numerical  stability  of  simple 
time  marching  schemes  used  in  elastic-viscoplastic  analysis.  The  Malvern 
and  Norton  models  were  studied  among  others.  For  the  Malvern  model  the 
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following  maximum  time  step  size  was  determined 


At  +  ±0+ v)  <r 

a 'm  -  3  £  (III-28) 

where  v  is  Poisson's  ratio  and  all  other  parameters  are  as  defined  earlier. 
For  the  Norton  model  the  maximum  time  step  was 


~  3  )!c£(3rta-‘ 


(III-29) 


In  general  dt^  is  several  orders  of  magnitude  larger  than  dt  •  If  the 
Malvern  model  or  the  Norton  model  are  used  separately  Cormeau  has  shown 
Eqs  III-28  and  III-29  to  work  well.  However,  if  the  Norton  and  Malvern 
models  are  coupled  together  for  a  more  complete  flow  law  a  method  is 
necessary  to  provide  the  maximum  time  step  during  transition  from  load 
up  (Malvern  dominated  phase)  to  creep  (Norton  dominated  phase).  In 
addition,  this  time  step  maximizing  scheme  is  all  the  more  required 
when  numerically  integrating  the  Bodner  equations  since  Cormeau's  analysis 
does  not  directly  apply  to  the  Bodner  viscoplastic  flow  law.  Consequently 
the  following  logic,  which  is  also  similar  to  the  MARC  program  [57],  was 
incorporated  into  VISCO. 

This  time  step  maximizing  logic  continually  tries  to  increase  the 
time  step  size  subject  to  a  stress  and  strain  constraint.  These  constraints 
are  the  allowable  amounts  of  change  in  stress,  a  and  strain,  Et0^> 
during  a  given  time  step  and  their  values  will  be  discussed  later  in  this 
section.  The  conputer  algorithm  is  based  around  parameters  and  P^ 
defined  as  follows 


(III-30) 


34 


(j<iy 

— : -  (III-31) 

^ to fa  I  *t.| 

\P1 - 2 - F* 

where  e  ^  =  Yc  +  e  +.50  Y  and  the  super  script  i  refers  to 

total  1  x  z  xy  r 

the  time  step.  Note  that  if  the  change  in  effective  stress  between  time 

steps  i  and  i-1  just  satisfies  the  stress  constraint  or  stress  tolerance, 

then  Eq  III-30  will  give  a  value  of  unity  for  P  .  Similarly,  if 

the  effective  plastic  strain  increment  for  time  step  i,  (dej^)*  ,  just 

satisfies  the  strain  constraint  or  strain  tolerance,  a  ,  ,  then  P  will 

tol  e 

equal  unity  from  Eq  III-31. 

The  parameters  and  P^  are  calculated  for  every  element  and  P  is 
set  equal  to  the  largest  one.  One  method  of  changing  the  time  step  size 
based  on  P  is 

dt‘  =  df'~l/p  (iix—32  ) 

Note  that  if  P  is  unity  no  change  in  the  time  step  size,  dt ,  occurs.  How¬ 
ever,  if  P  is  less  than  one,  dt*  is  greater  than  dt*  *  and  if  P  is  greater 
than  one  dt*  is  less  than  dt*  *.  In  the  case  of  P  being  greater  than  one, 
the  amount  of  change  in  stress  or  strain  has  exceeded  its  respective 
tolerance  and  recalculations  for  that  time  step  are  necessary  using  the 
reduced  time  step  size  from  Eq  III-32. 

To  avoid  several  successive  recalculation  steps  that  develop  when  P  is 
greater  than  one  and  Eq  III-32  is  used,  the  following  substitute  for  Eq 
III-32  was  employed. 
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el  t  ' 

-  .8  dt' 

/  p 

If 

P 

oft '  = 

-+• 

1 

Ti 

.8  ±  P 

/ 

ol  +  '  - 

I.Z  S’  ol+" 

1  If 

.£S  £  P 

< 

.  8 

eft'  ~ 

/.  5  clt 

If 

P  < 

.  b 

S' 

(Ill— 33) 


Note  that  Eq.  III-33  reduces  the  time  step  size  more  than  Eq .  III-32  if  P 
is  greater  than  one  but  if  P  is  less  than  one  Eq .  III-33  increases  the 
time  step  size  slower  than  does  Eq .  111-32.  Both  of  these  differences 
between  Eq.  III-33  and  Eq .  III-32  tend  to  reduce  the  number  of  calculation 
steps  required. 

Determination  of  the  values  for  the  stress  and  strain  tolerances, 

0^^  and  etol  respectively,  was  accomplished  by  employing  VISCO  to  analyze 
a  plate  with  a  V-notch.  This  particular  geometry  was  chosen  since  it  has 
a  high  stress  concentration  around  the  notch  and  is  therefore  somewhat 
similar  to  a  plate  with  a  crack  which  is  the  geometry  of  ultimate  interest 
in  this  research  effort.  However,  the  V-Notch  geometry  can  be  modeled  with 
far  less  elements  than  a  cracked  plate  requires  and  still  compare  with 
other  V-notch  solutions  in  the  literature  [40].  Figure  III-5  shows  the 
finite  element  mesh  employed  for  the  V-notch  plate  .  Only  one  quarter 
of  the  plate  is  modeled  due  to  symmetry.  Element  sizes  were  made  smallest 
near  the  V-notch  in  order  to  capture  the  stress  concentration  there. 

A  total  of  182  constant  strain  triangular  elements  were  employed  which 
is  somewhat  less  than  Yamada  et.  al.  [40]  used.  However,  good  agreement 
with  Yamada  was  achieved  and  this  will  be  demonstrated  later  in  the 
Validation  Examples  section. 
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Element  Number 


Figure  III-5.  V-Notched  Plate  Model 


MALVERN  MODEL.  An  elastic-perfectly  visco  plastic  (i.e.  o(c^)  = 

constant)  plane  stress  analysis  was  carried  out  using  the  Malvern  model. 

The  ratio  of  elastic  modulus  E  to  the  yield  stress  was  288  and  the  ratio 

of  applied  remote  stress  to  yield  stress  was  0.593.  Poisson's  ratio 

-4 

was  0.2.  The  load  was  applied  in  nondimensional  time  of  Jj^t  =  10 

where  Y  is  the  fluidity  constant  in  the  Malvern  model. Eq.  III-4. 

P 

Results  were  taken  after  all  stress  and  strain  rates  were  zero.  This 

may  be  defined  as  the  steady  state  condition  which  was  observed  to 

occur  after  t  =  0.4. 

P 

Several  analyses  of  the  V-notch  were  performed  using  different 

values  for  the  stress  and  strain  tolerances.  Examination  of  the  results 

found  that  between  displacements  near  the  notch,  plastic  strain  near 

the  notch,  and  total  plastic  strain  energy  in  the  model,  the  last  one 

was  most  sensitive  to  stress/strain  tolerance  variations.  Table  III-l 

displays  the  percent  variation  in  plastic  strain  energy  as  stress/strain 

tolerances  are  varied.  The  percent  variation  is  relative  to  the  plastic 

strain  energy  calculated  when  a  ,  =  .01  and  e  ,  =  .01. 

to 1  tol 

It  was  noted  that  the  stress/strain  tolerances  that  kept  the  percent 
variation  in  plastic  strain  energy  around  one  percent  also  kept  the  time 
step  size,  during  most  of  the  computing,  under  Cormeau's  critical  value, 
d  t^,  in  Eq.  III-2B.  The  amount  of  computer  time  increases  rapidly 
as  the  tolerances  are  reduced  to  a  value  of  .01.  A  good  compromise 
between  computer  time  required  and  apparent  accuracy  from  Table  III-l 
is  a  stress  tolerance  of  .03.  Note  that  for  this  stress  tolerance,  the 
strain  tolerance  can  be  relaxed  all  the  way  to  .20  with  little  change  in 
plastic  strain  energy.  These  results  are  in  agreement  with  the  recommendations 
for  stress/strain  tolerances  in  the  MARC  program  [57]. 
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BODNER  MODEL.  A  similar  suress/strain  tolerance  investigation  was 
performed  with  the  Bodner  material  model  in  plane  stress.  The  V-notch 
model  in  Fig  IIX-1  was  again  employed.  The  material  constants  used  for 
the  Bodner  equations  were  those  matched  to  IN-100  and  given  previously. 

The  ratio  of  applied  remote  stress  to  yield  stress  was  again  .593 
which  is  based  on  140  ksi  for  the  initial  yield  stress.  Maximum  load 
was  reached  in  10  seconds  and  plastic  strain  energy  results  were  then 
taken  similar  to  the  Malvern  model  investigation.  Again  several  analyses 
of  the  V-notch  were  performed  using  different  combinations  of  stress 
and  strain  tolerance  values.  The  results  for  the  Bodner  model  are  given 
in  Table  III-2.  Comparing  Table  III-l  and  III-2  shows  that  the  Bodner 
model  develops  less  plastic  strain  energy  variation  for  the  given  stress/ 
strain  tolerances  than  does  the  Malvern  model.  This  comparison  indicates 
that  higher  stress/strain  tolerances  could  then  be  used  for  the  Bodner 
model.  However,  if  the  analyses  are  continued  into  the  creep  phase, 
it  was  observed  that  the  stress  vali  es  tend  to  oscillate  somewhat  when 
they  should  be  monotonically  relaxir g  near  the  notch  tip.  This  oscillation 
was  fairly  well  damped  out  when  a  stress  tolerance  of  0.03  was  used. 

Once  this  small  value  for  the  stress  tolerance  is  chosen  the  strain 
tolerance  has  little  effect  as  long  as  its  greater  than  or  equal  to 
the  stress  tolerance. 

The  average  computation  time  for  the  Bodner  model  was  about  twice 
that  required  for  the  Malvern  ncdel.  The  average  central  processor 
time  required  on  the  CDC-6600  was  150  sec  for  this  200  degree  of  freedom 
problem  with  the  Bodner  model. 

VALIDATION  EXAMPLES 

The  following  five  examples  were  performed  to  demonstrate  the  validity 
of  the  VISCO  computer  program.  The  first  four  examples  employ  the  Malvurn 
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model  to  compare  with  published  time  independent  elastic-plastic  solutions. 
The  Bodner  model  is  then  tested  in  the  fifth  example  by  comparing 
its  behavior  with  the  results  from  coupling  the  Malvern  and  Norton  flow 
laws  together. 

INFINITE  SHEET  WITH  PRESSURIZED  HOLE.  The  first  example  employs  the 
Malvern  model  to  analyze  the  infinite  sheet  with  a  pressurized  hole. 
Comparisons  are  made  with  the  exact  elastic-plastic  strain  history 
independent  deformation  solution  by  Hsu  and  Forman  [70].  The  material 
properties  for  this  plane  stress  solution  were  E  =  10.5  x  10^  psi, 
v  =  0.5,  and  a  Ramberg-Osgood  stress-strain  equation 


€  = 


for 


<*1 


(III—  34 ) 


where  6  is  the  total  strain  (i.e.  elastic  plus  plastic),  oo  is  the 

initial  yield  stress  (aQ  =  55000  psi)  and  n  determines  the  shape  of  the 

curve  (n  =  9).  Eq .  III-34  must  be  inverted,  so  that  yield  stress  C  is  a 

function  of  plastic  strain  ,  in  order  to  be  used  in  the  Malvern  model. 

The  first  step  in  this  inversion  process  was  to  split  the  total  strain  c 

into  its  elastic,  ,  and  plastic,  eP  ,  parts  and  then  set  the  elastic 
E 

strain,  e  ,  equal  to  a/E.  Eq .  III-34  may  then  be  rewritten  as 


€  = 


<r_ 

E 


(III— 35) 


Solving  for  a  results  in 


crn  =  (V  +  E  £P) 


(II 1-36) 
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and  taking  the  n 


n+h 


root  of  Eq.  III-36 


cr  = 


(III-37) 


The  stress,  cr,  on  the  left  hand  side  of  Eq .  III-37  will  now  be  redefined 
as  the  strain  hardening  yield  stress  a  so  that,  approximately, 

—  HrJ.  /  _  .  tr  y  P 


cr  — 


C~  *0 


(«“  +  E  £p  )  "  -f.r  I  o- 1  -  <rc 


(111-38) 


Note  when  a  equals  oq  and  eP  is  zero,  Eq.  III-38  is  identically  satisfied. 
Also  when  o  is  greater  than  a  ,  the  plastic  strain  will  be  greater  than 
zero  and  the  product  E  e  Pwill  have  the  major  effect  on  o,  which  is  desired. 

The  finite  element  mesh  for  modeling  the  infinite  plate  was  similar 
to  Fig.  7b,  however,  132  triangular  elements  were  used.  The  outer  radius 
to  inner  radius  ratio  of  the  element  mesh  was  15  which  was  assumed  to 
approximate  an  infinite  radius  plate. 

Figure  III-6  shows  the  radial  and  tangential  stress  profiles  for 
three  internal  pressure  ratios,  P/oq.  Hsu  and  Forman  indicated  that, 
in  consideration  of  Budiansky's  criterion  for  the  acceptability  of 
deformation  theory,  their  infinite  plate  solution  should  not  disagree 
greatly  with  an  incremental  elastic-plastic  solution.  One  may  observe 
from  Fig.  III-6  a  close  approximation  between  the  techniques  thus 
lending  validity  to  the  incremental  plane  stress  approach  incorporated 
within  VISCO. 

THICK  CYLINDER  WITH  INTERNAL  PRESSURE.  The  second  example  employs 
the  Malvern  model  to  analyze  a  thick  cylinder  with  internal  pressure. 

In  this  case  comparison  is  made  with  a  non— finite  element  solution  by 
Hodge  and  White  [71]  who  again  used  deformation  theory.  This  was  an 
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RADIAL  &  TANGENTIAL  STRESSES/INITIAL  YIE^D  STRESS  £ 


P/«  =  1.3  , 
(tangential) 


VISCO  Results 


E  =  10.5E06  psl 
a;  =  55000.  psl 
~V  =  0.5 


Figure  III-6.  Radial  &  Tangential  Stresses  For  the  Infinite 
Plate  With  Pressurized  Hole 
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elastic-perfectly  plastic  plane  strain  analysis.  The  ratio  of  the  elastic 
modulus  to  yield  stress  was  E/a  =  190.9  with  v  =  0.33.  The  ratio  of  outer  radius 
to  inner  radius  for  the  cylinder  is  two.  Figure  III-7b  shows  the  finite  element 
mesh  used  to  model  a  symmetric  section  of  the  thick  cylinder.  The  mesh 
incorporates  40  triangular  elements.  Figure  III-7a  shows  the  tangential  stress 
profile  for  a  pressure  ratio,  P/a,  of  0.76.  Both  the  elastic  and  the  elastic- 
plastic  profiles  are  given.  There  is  good  agreement  with  Hodge  and  White  and 
thus  validity  is  again  given  to  the  incremental  plane  strain  portion  of  the 
VI SCO  program. 

V-NOTCHHD  PLATE  IN’  TENSION.  The  third  example  employs  the  Malvern  model 
to  analyze  a  V-notched  plate  in  tension.  Comparison  is  made  with  another 
finite  element  solution  by  Yamada,  et.  al.  [40].  A  time  independent  tangent 
modulus  approach  was  used  by  Yamada  et.  al .  for  their  elastic-plastic  analysis. 
This  was  an  elastic-perfectly  plastic  plane  stress  solution  where  the  ratio  of 
the  elastic  modulus  to  yield  stress  E/o  =  666.7  and  v  =  0.3.  The  finite  element 
mesh  used  herein  was  the  same  as  given  previously  in  Fig  III-5  where  only  one 
quarter  of  the  plate  is  modeled  due  to  symmetry.  Element  sizes  were  made 
smallest  near  the  V-notch  in  order  to  capture  the  stress  concentration  there. 

A  total  of  182  constant  strain  triangular  elements  were  employed  which  is  some¬ 
what  less  than  used  by  Yamada  et.  al. 

Figure  III-8  shows  two  y-component  of  stress  profiles  for  the  minimum 
section  of  the  V-notched  plate.  Note  for  the  applied  stress  of  a/a  =  0.185, 
which  was  Yamada' s  initial  yield  load,  the  present  results  from  VISC0  agree 
very  well  except  for  the  elements  immediately  at  the  notch.  The  high  elastic 
stress  concentration  near  the  notch  was  not  modelled  as  well  by  the  present  VIE  CO 
analysis  due  to  larger  element  sizes  being  used  in  the  notch  vicinity.  However, 
with  yielding  at  the  applied  stress  of  o/c  =  0.584,  comparison  with  Yamada's 
results  at  the  notch  are  much  better.  Note  that  plastic  action  diminishes  stress 
gradients  and  thus  fewer  elements  are  needed  for  an  elastic  solution  can  be  used 
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TANGENTIAL  STRESS 


Figure  1 1 1— 7 .  Thick  Cylinder  With  Internal  Pressure,  (a)  Tangential  Stress 
Profile,  (b)  Finite  Element  Mesh 
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&  VISCO  Elastic-Plastic  Results 


Figure  1 1 1  —  8 .  Stress  Profile  For  V-Notehed  Plate  In  Tension 
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Figure  III-9,  Plastic  Zone  In  V-Notched  Plate 
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to  produce  good  stress  profiles  after  yielding  occurs.  However,  it  should  be 
kept  in  mind  that  plastic  action  does  not  diminish  the  strain  gradient  like  the 
stress  gradient  and  therefore  strain  profiles,  even  after  yielding,  will  be 
quite  sensitive  to  element  sizing. 

Figure  III-9  displays  the  finite  element  mesh  for  the  V-notched  plate  with 
those  elements  left  out  that  have  exceeded  the  yield  stress  for  the  applied 
stress  of  a/a  =  0.584.  The  absence  of  these  elements  thus  describes  the  plastic 
zone  size  compares  well  with  that  of  Yamada.  Numbers  within  the  elements  of 
Fig  III  -9  reflect  each  elements  effective  stress  as  a  percentage  of  yield 
stress . 


CRACKED  THREE  POINT  BEND  SPECIMEN.  The  fourth  example  employs  the  Malvern 
model  to  analyze  a  cracked  three  point  bend  specimen  as  shown  in  Fig  III-10. 
This  will  be  the  same  problem  as  was  solved  with  ten  different  computer  codes 
for  an  ASTM  analytical  round  robin  and  documented  by  Wilson  and  Osias  [72]. 

This  was  a  plane  strain  crack  problem  where  E=  31-59  x  10^  psi,  v=  0.30,  and 
the  following  Rambe rg-Osgood  equation  was  used  for  the  stress-strain  curve 


where  e  is  the  total  strain  (i.e.  elastic  and  plastic), 
is  120  x  10^  psi.  To  make  Eq  III-39  compatible  with  the 
stress  must  be  written  as  a  function  of  plastic  strain, 
tracting  the  elastic  strain  (also  equal  to  a/E)  from  both 

solving  for  the  stress  a,  and  redefining  a  as  a 

_L 
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n  is  10,  and  BQ 
Mai vern  model , 
Therefore,  after  sub¬ 
sides  of  Eq  HI-39, 
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Due  to  symmetry  only  one  half  of  the  three  point  bend  specimen  was 
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Figure  III-10.  Three  Point  Bend  Specimen  With  Crack 


0 


5 


modeled  with  the  finite  element  mesh  in  Fig.  III-ll.  This  particular 

pattern  for  the  element  mesh  was  used  by  Ohtani  and  Nakamura  [61]  for  a 

center  cracked  plate.  Note  how  the  element  sizes  are  reduced  as  indicated 

in  Fig.  III-ll  by  arrows  to  the  first  and  second  reduction.  This  pattern 

provides  for  an  unlimited  number  of  element  size  reductions  while  also 

maintaining  good  element  aspect  ratios  (e.g.  from  1  to  0.5)  and  insuring 

no  two  neighboring  elements  differ  in  size  by  more  than  a  factor  of  2. 

Accordingly,  each  reduction  cuts  the  preceding  element  size  in  half. 

Eight  element  size  reductions,  incorporated  in  the  Fig.  III-ll  mesh, 

step-wise  reduced  the  0.20  in.  elements  at  the  upper  boundary  to  a 
-4 

7.8125  x  10  in.  element  at  the  crack  tip.  This  crack  tip  element  size 
was  slightly  smaller  than  the  smallest  used  in  reference  [72].  The  total 
number  of  elements  in  Fig.  III-ll  was  584  with  388  nodes.  Figures  III-12, 
III-13,  and  III-14  show  the  results  of  VISCO  compared  to  10  time  independent 
elastic-plastic  solutions  documented  by  Wilson  and  Osias.  Eight  of 
these  solutions  fall  within  lines  II  and  III  in  Figs.  III-13  and  III-14. 

The  other  two  solutions  fall  near  line  I  in  these  two  figures.  In  all 
cases  the  present  results  from  VISCO  fall  within  the  ASTM  round  robin 
range  as  shown  in  these  figures.  Therefore  the  element  sizing  and 
arrangement  in  Fig.  III-ll  as  used  herein  with  VISCO  appears  to  give 
good  results  for  a  complex  problem  that  includes  a  crack  and  loads  that 
induce  bending. 

CENTER  CRACKED  PLATE  WITH  BODNER  MODEL.  The  fifth  example  employs 
the  Bodner  model  to  analyze  a  center  cracked  plate  ,  Comparison  will  be 
made  with  results  from  a  similar  analysis  using  the  Malvern  model  coupled 
with  Norton's  Law  as  follows 


-<rK) 

L*C  4^  .■*  «cfr 
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Figure  rri-11.  Three  Point  Bend  Finite  Element  Model 


Carck  Mouth  Opening  Displacement  vs.  Load  For  Three  Point  Model 


This  Malvern-Nor ton  combination  is  a  superposition  approach  suggested 
by  Zienkiewicz  [10]  to  model  in  a  unified  sense  both  initial  load  up 
viscoplasticity  and  creep  under  sustained  load.  Also,  to  assess  the 
contribution  of  pure  secondary  creep  in  the  Bodner  model,  comparison 
will  be  made  to  an  analysis  using  only  Norton’s  Law  for  the  viscoplastic 
material  model. 

The  material  properties  will  be  those  matched  to  IN-100  at  1350°F. 

The  constants  for  the  Bodner  model  will  be  those  previously  given  for 
IN-100.  The  yield  stress  as  a  function  of  plastic  strain,  o(e^),  needed 
for  the  Malvern  model  will  be  a  multilinear  fit  to  the  experimental 
stress-strain  curve  in  Fig.  III-2  and  given  as  follows 

f  I  3  o.  (l.  +  7(> .  2  )  K  s  i  for  6% <  .  ooZZZ 

)  (HI-42) 

a~  ((*  )  -  A  /S2.4  for  .oozzz*  Ce  <  •<>  i 

(  4  KS1  f°r 


The  fluidity  constant  y^  for  the  Malvern  model  will  be  given  a  value 

of  ,58  sec  \  This  value  for  y^  is  sufficiently  high  such  that  the 

given  a  function  is  followed  very  closely  (e.g.  recall  Fig.  III-l  for 

y  =  °°)  .  The  constants  for  Norton's  Law  will  be  y  =  3.7394  x  10 
P  c 

(psi)  ^‘^/sec  and  8  *  10.64.  These  values  were  determined  by  matching 
the  initial  IN-100  creep  behavior  in  Figs.  III-3  and  III-4  as  discussed 
in  the  Norton's  Law  for  secondary  creep  section. 

Due  to  symmetry  only  one  fourth  of  the  center  cracked  p]ate  is 


modeled  by  the  finite  element  mesh  given  in  Fig.  III-15.  The  plate’s 
height  is  2.8  in.,  width  is  1.  in.,  and  thickness  is  0.3  in.  The  crack 
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Figure  1 1 1  —  1 5 .  Center  Cracked  Plate  Finite  Element  Model 
a/W  =  0.1367 
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length,  2a,  is  .2734  in.  or  a/W  equal  to  .1367.  The  numbers  inside  the 
elements  indicate  a  total  of  355  elements  were  employed.  The  total  number 

of  nodes  is  211.  The  triangular  elements  around  the  crack  tip  have  a 
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height  and  base  dimension  of  7.8125  x  10  in.  Further  discussion  of 
this  element  mesh  will  be  presented  in  the  Applications  section  since 
this  particular  mesh  was  also  employed  there  to  simulate  the  experimental 
program.  A  maximum  stress  of  36320  psi  will  be  applied  to  the  upper 
boundary  of  the  center  cracked  plate  in  5  seconds  and  then  held  constant 
for  an  elapsed  time  of  1000  seconds.  This  applied  stress  level  corresponds 
to  a  load  level  also  used  in  the  experimental  program  to  be  simulated  later. 

A  comparison  of  the  behavior  of  these  three  material  models  is 
given  in  Figs.  III-16,  III-17,  and  III-18.  The  effective  stress  and  plastic 
strain  in  Figs.  III-16  and  III-17  are  from  the  element  at  the  crack  tip 
which  had  the  highest  elastic  stress  concentration  factor. 

The  stress-strain  behavior  in  Fig.  III-16  occurred  over  a  time  period 
of  1000  sec.  The  values  to  the  left  of  6%  strain  occurred  in  approximately 
5  sec  (the  load  up  period)  whereas  to  the  right  of  6%  strain,  stress 
relaxation  and  redistribution  is  taking  place  over  approximately  1000  sec 
of  sustained  1-^ad  creep  behavior.  Note  how  the  Norton  model  behaves 
during  the  load  up  period  (i.e.  strains  less  than  6%).  Due  to  its  slow 
creep  response,  the  effective  stress  from  the  Norton  model  approaches 
the  values  that  would  occur  during  load  up  in  an  elastic  analysis.  However, 
the  Norton  model  relaxes  the  stress  down  to  values  very  similar  to  the 


other  material  model  stress  values.  In  general  the  Bodner  model  behaves 
very  similar  to  the  Malvern-Nor ton  combination  both  during  load  up  and 
the  sustained  load  creep  periods. 
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Figure  III-16.  Stress  vs.  Strain  At  The  Cracktip  For  Bodner,  Malvern- 
Norton,  and  Norton  Material  Models 
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Figure  1 1 1— 1 7 .  Effective  Stress  &  Strain  vs.  Time  At  The  Cracktip  For 
Bodner,  Malvern-Norton,  and  Norton  Material  Models 
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Figure  III-17  shows  the  time  dependent  behavior  of  the  effective 
stress  and  plastic  strain  from  the  same  crack  tip  element.  After  approximately 
200  sec.,  all  three  material  models  develop  nearly  indentical  effective 
stress  values.  Plastic  strain  behavior  with  time  is  also  very  similar 
for  the  three  models  except  for  under  100  sec.  of  time.  The  difference 
in  percent  plastic  strain  between  the  three  curves  developed  primarily 
from  different  plastic  strain  rates  during  the  load  up  phase  and  remains 
fairly  constant  for  time  greater  than  200  sec. 

Figure  III-18  shows  how  the  crack  mouth  displacement  increases 
with  time  after  the  maximum  load  is  achieved.  The  location  of  the  crack 
mouth  displacement  is  indicated  in  Fig.  III-18  to  be  0.050  in.  from  the 
vertical  centerline.  Again  the  three  curves  are  very  similar  after  200  sec. 
and  their  separation  is  due  primarily  to  dissimilar  behavior  for  time 
under  200  sec.  The  Norton  model  displays  more  displacement  after  maximum 
load  than  the  other  models  since  it  is  effectively  making  up  for  its 
slower  plastic  strain  rates  and  associated  displacements  during  the  load 
up  phase.  This  apparently  is  also  true  when  comparing  the  Malvern-Nor ton 
to  the  Bodner  model  displacement  curve,  however,  here  the  difference 
is  much  less. 

Therefore,  based  on  these  comparisons  with  somewhat  similar  material 
models  used  to  analvze  the  center  cracked  plate,  the  Bodner  model  within 
VISCC  is  considered  to  be  working  well.  For  further  discussion  and 
comparison  of  these  material  models  sl^  reference  [77]. 

SUMMARY  OF  FINDINGS 

The  following  findings  were  determined  from  the  preceding  validation 
examples  and  the  general  development  of  the  VISC0  program. 

1.  The  stress  tolerance,  a  . ,  which  controls  the  variable  time 

tol 

step  size  in  VISCO's  numerical  time  integration  algorithm,  provides 
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1 

I 

good  results  for  reasonable  computer  time  requirements  when  set  \ 

i 

to  the  value  of  0.03.  The  strain  tolerance,  e  ,  ,  has  little  * 

tol 

effect  as  long  as  its  greater  than  or  equal  to  the  stress  tolerance. 

This  finding  pertains  to  the  Bodner,  Malvern,  and  Norton  material 
models  . 

2.  The  VISCO  program,  while  employing  the  Malvern  material  model 
agrees  well  with  so-called  exact  elastic-plastic  deformation 
solutions  both  in  plane  stress  and  plane  strain  conditions. 

Agreement  is  also  good  with  time  independent  elastic-plastic  ; 

*i 

finite  element  solutions  in  a  V-notched  plate  and  a  cracked 

three  point  bend  specimen  whose  results  came  from  an  ASTM  -f 

analytical  round  robin. 

3.  The  finite  element  mesh  pattern  in  Fig.  III-ll  works  well  for 

modeling  cracked  plates.  This  pattern  conveniently  provides 

for  an  unlimited  number  of  element  size  reductions  to  capture 

the  crack  tip  singularity  while  also  maintaining  good  element 

aspect  ratios  and  minimizing  the  total  number  of  elements  required. 

-4 

4.  A  crack  tip  element  size  to  specimen  width  ratio  of  7.8125  x  10 
in  a  cracked  three  point  bend  specimen  provided  good  agreement 
with  an  ASTM  analytical  round  robin. 

5.  The  Bodner  material  model  has  been  shown  to  behave  similarly  to 


the  Malvern-Nor Lon  model  for  both  load  up  and  sustained  load  creep 
stages.  Also,  for  time  greater  than  200  sec.  after  load  is  applied, 
the  pure  secondary  creep  Norton  law  behaves  similar  to  the 
Bodner  model.  In  general,  these  similarities  in  material  model 
behavior  should  be  true  for  most  any  metal  however,  the  indicated 
200  sec.  time  delay  between  the  Norton  and  Bodner  models  pertains 
specifically  to  the  IN-100  alloy  at  1350°F. 


i 


These  findings  therefore,  support  the  validation  of  the  VISCO 
computer  program  arid  provide  some  of  the  required  details  for  applying 
the  VISCO  program  with  the  Bodner  model  to  creep  crack  growth  simulation. 
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IV  HYBRID  EXPERIMENTAL- NUMERICAL  PROCEDURE  TO 
ANALYZE  CREEP  CRACK  GROWTH 

The  simultaneous  use  of  experimental  data  from  crack  growth  tests 
and  a  theoretical  model  of  the  experimental  cracked  specimen  has  been 
termed  the  Hybrid  Experimental-Numerical  procedure  by  Kobayashi  [16]. 

One  example  of  this  procedure  would  be  to  grow  a  crack  at  experimentally 
determined  rates  through  the  theoretical  model  of  the  experimental  specimen. 
Then,  from  the  results,  one  could  seek  out  potential  crack  growth  criteria 
which  hopefully,  during  crack  growth,  display  themselves  as  fairly 
constant  parameters  independent  of  both  crack  length  and  load  intensity 
(e.g.  stress  or  strain  at  the  crack  tip,  crack-opening  displacement,  etc.). 

In  the  present  investigation,  the  Hybrid  Experimental-Numerical 
procedure  required  the  theoretical  model  to  track  experimental  displacement 
rates  rather  than  crack  growth  rates.  The  theoretical  model  consisted 
of  the  VISCO  program  employing  the  Bodner  constitutive  equations  and  a 
finite  element  mesh  of  the  experimental  specimen.  The  experimental  dis¬ 
placements  were  either  measured  near  the  crack  tip  or  near  the  vertical 
centerline  (crack  mouth)  of  the  center  cracked  plate  test  specimens  as  shown 
in  Fig  IV-1.  The  displacements  were  measured  continuously  with  time 
by  an  optical  interferometric  displacement  measurement  technique  developed 
by  Sharpe  [74],  Figure  IV-2  shows  a  typical  experimental  displacement 
versus  time  curve.  The  present  numerical  procedure  required  the  crack 
to  grow  through  the  theoretical  model  such  that  the  displacements  due  to 
elastic-plastic  behavior  and  crack  growth  added  up  to  the  experimental 
displacements  as  time  progressed.  Therefore,  crack  length  versus  time 
became  a  product  of  the  present  analysis  rather  than  an  input. 

CRACK  LENGTH  VERSUS  TIME 

Ideally  for  this  investigation,  the  experimental  data  should  be  in 
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Figure  IV-2.  Displacement  Under  Constant  Load 


the  form  of  crack  length  versus  time  rather  than  displacement  rates. 
Unfortunately,  it  was  nearly  impossible,  experimentally,  to  measure  crack 
length  as  a  function  of  time  with  any  degree  of  precision  or  reliability  on 
the  surface  since  the  total  amount  of  crack  growth  in  these  tests  was  extremel 
small  (e.g.  100  microns)  and  since  the  creep  crack  would  grow  internally 
or  tunnel  without  associated  surface  crack  growth. 

Another  attempt  to  measure  crack  growth  rate  indirectly  employed 
the  elastic  compliance  method  first  demonstrated  by  Clarke  [75].  This 
method  utilizes  the  change  in  elastic  compliance  of  the  specimen  with 
time  and  then  with  the  aid  of  a  compliance  versus  crack  length  relationship 
based  on  linear  elastic  fracture  mechanics,  crack  growth  with  time  may  be 
determined.  Figure  IV-3  shows  schematically  an  experimental  creep  crack 
growth  load  versus  time  history  designed  to  provide  discrete  values  of 
compliance  at  selected  time  intevals.  The  load  is  reduced  approximately 
20%  and  then  restored  as  shown  at  times  t  through  t^  to  provide  load  dis¬ 
placement  data  at  various  times  during  the  test.  Figure  IV-4  shows  a 
typical  set  of  load  displacement  data  for  increasing  times  t^  to  t^  in  a 
creep  crack  growth  test.  Compliance  is  defined  as  displacement  divided 
by  load  and  therefore  the  slope  of  each  line  in  Fig  IV-4  represents  the 
compliance  at  each  particular  time.  Note  that  compliance  is  shown  to 

decrease  in  going  from  time  t.  to  t„  and  then  increase  from  time  t„  to  t,. 

i  1  l  4 

Comparing  this  behavior  to  a  typical  elastic  compliance  versus  crack  length 
curve  given  in  Fig  IV-5,  the  mathematical  implication  is  that  the  crack 
shortens  with  time.  Although  this  compliance  decrease/increase  behavior 
has  also  been  observed  by  Donat  [76],  no  known  experimental  data  supports 
any  physical  shortening  or  healing  of  the  crack.  However,  the  important 
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Figure  IV-3.  Experimental  Load  vs.  Time 


Figure  IV-4 .  Experimental  Compliance  Changes  With  Time 
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Figure  IV-5.  Elastic  Compliance  vs.  Crack  Length 
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implication  of  this  compliance  behavior  is  that  a  one  to  one  relation 
between  crack  length  and  compliance  does  not  exist  during  creep  crack 
growth.  Therefore,  the  elastic  compliance  method  cannot  be  used  directly 
to  measure  crack  length  with  time  in  creep  crack  growth  tests. 

MATCHING  EXPERIMENTAL  DISPLACEMENT  DATA 

In  each  Hybrid  Experimental-Numerical  application  of  the  VISCO 
finite  element  model  to  simulate  creep  crack  growth  tests,  the  y-dis- 
placement  of  the  node  which  was  closest  to  the  point  where  the  experimental 
measurement  was  made  was  monitored  with  time.  If  the  node's  displacement 
became  less  than  experimental  displacement  at  a  given  time  a  crack  tip 
node  would  begin  to  be  released  to  simulate  crack  growth  through  the  model. 
CRACK  TIP  NODE  RELEASE  METHODS 

Crack  tip  nodes  were  released  in  one  of  two  different  methods.  The 
first  method  releases  the  node  and  totally  unloads  it  in  5  sec.  The 
second  method  unloads  the  current  crack  tip  node  linearly  with  time  over 
the  total  time  span  between  the  time  when  the  current  crack  tip  node  is 
begun  to  be  released  and  the  time  when  the  next  crack  tip  node  will  be 
released.  The  5  sec.  node  release  method  for  crack  growth  must  be  used 
when  matching  experimental  displacements  or  when  a  crack  growth  criterion 
is  used.  In  both  of  these  cases  when  certain  conditions  are  satisfied, 
the  crack  must  grow  so  a  node  is  released.  However,  when  the  current  node 
is  being  unloaded  it  is  not  known  when  the  next  crack  tip  node  will  be 
released  and  therefore  the  continuous  unload  method  cannot  be  used.  The 
5  sec.  unload  time  for  the  first  method  was  based  on  the  size  of  the  crack 
tip  elements  and  the  maximum  crack  growth  rates  occurring  in  the  creep 
crack  growth  test  data  (i.e.  maximum  crack  growth  rate  equals  element 
size  divided  by  5  sec).  If  a  crack  growth  rate  criterion  were  used  then 
it  could  be  determined  by  extrapolation  when  the  next  crack  tip  will  be 


released  and  thus  the  current  crack  tip  node  could  be  unloaded  in  a 


continuous  fashion  by  the  second  method.  The  second  node  release  method 
can  also  be  employed  when  all  node  release  times  are  specified  at  the 
beginning  of  the  computer  run  (e.g.  release  times  based  on  results  from 
a  prior  computer  run  using  release  method  one  and  matching  experimental  data)  . 

In  both  node  release  methods,  the  node  force  required  to  hold  the 
crack  tip  node  at  zero  displacement  would  be  calculated  from  the  stresses 
in  the  adjacent  elements  as  follows 

~  Z  dV-1) 

vo  I 

E| e »»i  e  v>  4s 

which  is  consistant  with  the  formulation  of  the  stiffness  matrix  in 
Appendix  A.  The  crack  tip  node  restraint  force,  f  ,  is  then  the  component 
of  { f }  perpendicular  to  the  crackline.  The  boundary  condition  on  the  node 
is  converted  from  zero  displacement  to  a  force  equal  to  f  .  This  force 

y 

f  is  then  removed  depending  upon  which  node  release  method  is  chosen, 
see  Fig  IV-6. 

The  change  of  the  crack  tip  node's  boundary  condition  from  dis¬ 
placement  to  a  force  boundary  condition  is  handled  very  conveniently  with 
the  Gauss-Seidel  iterative  linear  equation  solver  as  discussed  in  Appendix  B. 
Whenever  a  node  is  fixed  in  a  certain  direction,  its  equilibrium  equation 
in  that  direction  is  skipped  over  during  the  iterative  solution  procedure 
and  when  the  node  is  released  its  equilibrium  equation  is  included  within 
the  iterative  procedure.  No  lengthy  refactoriza  ion  of  the  stiffness  matrix 
is  required  for  these  boundary  condition  changes. 
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V  APPLICATION  OF  THE  HYBRID  EXPERIMENTAL 
NUMERICAL  PROCEDURE  TO  CREEP  CRACK  GROWTH 
The  first  objective  in  the  present  research  was  to  develop  a  method  for 
getting  crack  growth  behavior  solely  from  displacement  measurements  made  on  a 
cracked  specimen  under  constant  load  and  at  elevated  temperature.  This  objec¬ 
tive  is  extremely  important  since  small  but  significant  amounts  of  crack  growth 
can  not  otherwise  be  resolved  by  conventional  experimental  crack  measuring 
techniques.  The  second  objective  was  to  seek  out  crack  growth  criteria  based 
on  the  crack  growth  behavior  identified  from  the  present  work  on  the  first 
objective  and  by  examining  various  parameters  around  the  cracktip  in  the  the¬ 
oretical  model  (e.g.  stress,  strain,  crack  opening  displacement,  etc.). 

This  section  presents  the  results  of  applying  the  hybrid  experimental 
numerical  (HEN)  procedure  to  creep  crack  growth  in  IN-100  at  1350°  F.  The 
experimental  portion  of  the  procedure  consisted  of  displacement  versus  load 
(i.e.  compliance)  and  displacement  versus  time  test  data  reported  by  Sharpe 
[  15],  The  numerical  portion  of  the  HEN  procedure  consisted  of  the  VISCO 
finite  element  program  employing  the  Bodner  material  model.  The  material  con¬ 
stants  for  the  Bodner  model  were  those  given  in  Section  III. 

The  machining  specifications  for  the  specimens  used  in  the  experimental 
program  are  shown  in  Fig  V-l.  Only  the  center  uniform  cross  section  part  of 
the  specimen  was  considered  in  the  VISCO  simulation  and  due  to  symmetry  only 
one  quadrant  of  this  section  was  represented  by  the  finite  element  meshes 
given  in  Figs  V-2,  V-3 ,  and  V-4.  Each  of  these  meshes  represent  the  center 
cracked  plate  test  specimen  with  different  crack  lengths.  The  convergence  of 
these  meshes  has  been  verified  through  the  work  in  Section  III  and  further 
discussion  of  their  accuracy  will  be  provided  in  subsequent  paragraphs.  Fig- 
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res  V-2,  V-3,  and  V-4  have  half  crack  lengths,  a,  or  0.137  in.,  0.237  in., 
0.312  in.  respectively.  Due  to  the  variations  of  the  initial  crack  length 
through  the  thickness,  the  surface  measurement  cannot  be  used  directly  and 
thus  effective  initial  crack  lengths  must  be  determined.  These  effective 
initial  crack  lengths  were  determined  by  matching  experimental  load-dis¬ 
placement  data  along  the  crack  such  that  the  finite  element  model  displayed 
elastically  the  same  compliance  as  the  experimental  specimens  did.  For 
efficiency  sake  it  has  been  found  realistic  to  use  the  same  element  mesh  for 
slightly  different  initial  crack  lengths. 

The  VISCO  finite  element  program  has  both  plane  stress  and  plane  strain 
analysis  capability.  The  current  investigation  chose  plane  stress  as  re¬ 
ported  in  reference  [78]  where  theoretical  plane  stress  J  integral  values 
agreed  best  with  test  data  for  an  even  thicker  (  1  in.  thick  )  compact  ten¬ 
sion  specimen.  A  theoretical  model  must  display  realistic  compliance  be¬ 
havior  in  order  to  calculate  J  values  that  agree  with  test  data.  Likewise,  in 
the  present  research  realistic  compliance  behavior  is  a  necessity. 

Figures  V-2,  V-3,  and  V-4  display  the  same  general  pattern  of  elements 

which  worked  well  for  the  three  point  bend  specimen  in  Section  III.  The  ele- 

-4 

ments  at  the  cracktip  have  a  height  and  width  of  7.8125  x  10  in.  This  size 
cracktip  element  in  combination  with  the  given  general  element  pattern,  pro¬ 
vided  for  355,  278,  and  362  elements  respectively  in  the  three  figures. 

Figures  V-2  and  V-4  provided  20  uniform  elements  ahead  of  the  cracktip  for 
subsequent  crack  growth  whereas  Fig  V-3  had  eight  uniform  elements.  Figure 
V-5  shows  the  expanded  element  layout  around  the  cracktip  from  Fig  V-3.  This 
region  of  uniform  elements  ahead  of  the  initial  cracktip  avoids  unrealistic 
changes  in  compliance,  as  the  crack  grows  through  the  model,  that  can  develop 

when  nonuniform  element  sizes  are  used.  The  number  of  uniform  elements 
incorporated  ahead  of  the  cracktip  is  a  compromise  between  anticipated  crack 
growth  and  computer  time  required. 
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The  elastic  compliance  of  the  specimens  modelled  by  these  three  meshes 
in  VISCO  was  compared  to  an  empirical  solution  by  Eftis  and  Liebowitz  [79] 


in  the  form  of 


CSC 


c:osh“‘[s*cC-57')'] 


(v-i) 


This  comparison  is  shown  in  Fig  V-6,  with  the  symbols  a,v,and  W  as  defined 
in  the  figure.  Also  shown  is  nondimens ionali zed  compliance  versus  crack 
length  for  locations  off  the  vertical  centerline  and  relatively  near  the 
cracktip.  These  locations  correspond  to  the  displacement  measurement  loca¬ 
tions  used  in  the  experimental  program.  It  can  be  seen  that  the  centerline 
compliance  agrees  very  well  with  published  results.  Also  note  that  if  the 
off  centerline  results  are  linearly  extrapolated  (i.e.  dashed  line)  back  to 
the  Effis  and  Liebowitz  curve  they  intersect  at  a/W  values  which  correspond 
to  their  distance  behind  the  cracktip.  These  results  all  support  the  validity 
of  the  finite  element  meshes  used  in  the  present  investigation. 

Figure  V-7  shows  elastic  crack  opening  displacement  profiles  using  VISCO 
and  the  mesh  in  Figure  V-2  compared  to  the  Westergaard  equation  for  elastic 
displacements  around  the  cracktip.  The  following  form  of  the  Westergaard 
equation  is  restricted  to  plane  stress  displacements  behind  the  cracktip  and 
on  the  crack  surface  (i.e.  crack  opening  displacement) 


where  r  is  the  distance  behind  the  cracktip  and  is  the  mode  I  elastic  stress 
intensity  factor  which  for  the  center  cracked  plate  can  be  written  as  [73] 
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Figure  V-6.  Comparison  of  Compliance  From  VISCO  With  Eftis  & 
Llebowitz  For  The  Center  Cracked  Plate 
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where  again  a  and  W  are  che  half  crack  length  and  specimen  width  respectively. 
Mote  that  agreement  with  the  Westergaard  equation  is  quite  good  here  for  the 
results  from  VISCO  using  the  mesh  with  the  shortest  crack  length.  Fig  V-2. 

The  same  cracktip  element  size  was  used  in  all  three  crack  length  models. 
Therefore  agreement  with  Eq  V-2  would  even  be  better  for  the  longer  crack 
length  models  since  the  ratio  of  element  size  to  crack  length  will  be  smaller. 
It  has  been  shown  through  elastic  finite  element  convergence  studies  that  as 
this  ratio  of  cracktip  element  size  to  crack  length  decreases,  accuracy  in¬ 
creases  [81]. 

Another  consideration  supporting  this  size  cracktip  element  (i.e.  7.8  x 
-4 

10  in.)  is  that  IN-100  grains  are  approximately  the  same  dimension  [13]. 

It  may  be  argued  that  the  finite  element  method  is  a  continuum  analysis  tool 
and  that  accuracy  should  only  improve  as  elements  are  refined.  However, 
realistically  a  continuum  does  not  exist  at  and  below  the  grain  size,  espe¬ 
cially  around  the  cracktip.  Thus  incorporation  of  elements  smaller  than  the 
grain  size  might  be  unrealistic.  Furthermore,  it  is  postulated  that  the 
physical  blunting  of  the  cracktip  that  can  be  associated  with  noncontinuous 
grain  structured  material  might  be  more  effectively  considered  by  the  finite 
element  model  used  herein  since  the  cracktip  elements  have  a  dimension  on  the 
order  of  a  grain  size. 

CRACK  GROWTH  PREDICTIONS 

A  summary  of  the  creep  crack  growth  test  data  reported  by  Sharpe  [15] 
la  given  in  Tables  V-l  and  V-2.  In  general  the  experimental  program  objectives 
were  to  generate  creep  crack  growth  data  in  the  center  cracked  plate  speci- 
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men  at  1350°  F  for  several  different  loads  and  cracked  lengths.  The  loads 
were  specified  in  terms  of  a  range  of  stress  intensity  factor  values  from 
approximately  15.0  to  35.0  ksi  V  in’  for  each  respective  initial  crack  length. 
With  these  objectives  in  mind  and  the  limited  number  of  test  specimens, only 
one  test  was  done  at  each  of  the  test  conditions.  In  the  process  of  de¬ 
veloping  the  experimental  procedure,  several  tests  did  not  result  in  good 
data  and  consequently  were  not  used  in  the  present  investigation  as  implied 
by  the  discontinuous  test  numbers  in  Table  V-l. 

The  initial  surface  crack  lengths,  ag,  were  measured  on  the  surface  of 
the  test  specimens  after  fatigue  precracking  and  prior  to  the  load  applica¬ 
tion  for  the  creep  crack  growth  tests.  Displacements  across  the  crack  were 
measured  at  the  indent  locations,  qualitatively  referred  to  as  mouth  or  tip 
locations  described  in  Section  IV.  In  each  of  these  creep  crack  growth  tests, 
the  load  was  applied  in  5  sec.  and  held  constant  for  the  given  test  duration 
time  (excluding  small  unload/reload  cycles  for  compliance  referred  to  in 
Section  IV) . 

Figure  V-8  is  a  photograph  of  the  fracture  surface  of  specimen  number 
2.  This  photograph  was  taken  after  breaking  open  the  center  cracked  plate 
specimen  following  creep  crack  growth  tests  5  and  6.  The  two  dark  bands  in¬ 
dicated  with  arrows  are  the  creep  crack  extensions  from  tests  5  and  6  where¬ 
as  the  remaining  fracture  surface  is  either  pre-test  or  post-test  fatigue 
crack  growth. 

Post-test  crack  length  measurements  were  made  on  the  fracture  surface 
as  described  in  Fig  V-9.  The  average  of  the  four  crack  length  measurements 
was  defined  as  the  experimental  crack  length  a^.  This  averaging  was  done  to 
smooth  out  crack  length  differences  due  to  asymmetric  crack  growth  and  vari- 


ations  through  the  thickness.  Table  V-2  gives  a  tabulation  of  the  initial 
crack  length,  ag,  and  am-  In  addition,  a  crack  length  determined  by  com¬ 
pliance,  a^,  and  the  initial  crack  length  used  in  the  finite  element  model, 
aQ,  is  given.  The  crack  length  was  determined  as  discussed  earlier  by 
varying  the  crack  length  in  the  finite  element  model  until  the  model's  com¬ 
pliance  matched  experimental  compliance  data. 

For  convenience  a  modified  form  of  the  above  compliance  technique  for 
crack  length  determination  was  also  employed  which  reduced  the  number  of 
element  meshes  for  different  crack  lengths  required.  This  modification  to 

the  compliance  method  made  use  of  the  experimental  compliance,  C  ,  and  the 

h 

finite  element  model  compliance,  C_,_,  for  a  crack  length  near  the  test 

r  h 

value.  To  determine  the  test  specimen's  crack  length,  the  difference  be¬ 
tween  CE  and  CpE  was  divided  by  the  rate  of  change  of  compliance  with  crack 
length  /Sc/ IS a  (i.e.  slope  of  curve  in  Fig  V-6)  as  follows  and  added  to  the 
model's  crack  length 


ac  -  Qo  + 


Cb  ~  CFj 
ACAa 


(V-4) 


where  c/  a  is  determined  from  the  model  for  nearby  crack  lengths.  The  com¬ 


pliance  CpE  pertains  to  a  model  crack  length  of  aQ. 

The  effective  stress  intensity  factor,  in  Table  V-2  was  calculated 

from  Eq  V-3  using  the  crack  length,  a^,  and  the  load  given  in  Table  V-l.  The 
stress  intensity  factor,  KpE>  is  also  calculated  for  convenience  from  Eq  V-3 
but  using  the  crack  length  aQ.  For  another  check  on  the  element  mesh,  the 
elastic  stress  intensity  factor,  K^,  was  also  calculated  based  on  J  integral 
values  determined  from  VISCO  for  several  paths.  For  linear  elastic  plane 


stress  behavior  [73] 


Kx  =  Al~ 


(V-5) 
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(Appendix  C  describes  the  VISCO  routine  for  calculating  the  J  integral). 

Figure  V-lOa  shows  a  scaled  drawing  giving  four  different  J  integral 
paths  used  in  VISCO  for  the  center  cracked  plate  specimen.  Figure  V-lOb 
shows  normalized  stress  intensity  factors  calculated  from  Eq  V-4  and  the  J 
values  along  paths  one  through  four  in  Fig  V-lOa .  These  J  values  were  from 
a  linear  elastic  VISCO  analysis.  For  linear  elastic  material,  the  J  in¬ 
tegral  is  theoretically  path  independent  and  this  path  independence  is 
demonstrated  in  Fig  V-10.  It  should  be  noted  that  the  good  agreement  be¬ 
tween  VISCO  results  and  Eq  V-3  in  Fig  V-lOb  also  indicates  that  the 
finite  element  mesh  (i.e.  Fig  V-2)  accurately  represents  the  center 
cracked  plate. 

Table  V-3  summarizes  the  basic  details  of  the  VISCO  computer  runs  which 
employed  the  HEN  procedure.  The  computer  run  numbers  designated  by  SI 
through  S7  will  frequently  be  referred  to  in  the  following  discussions.  As 
described  in  section  IV,  each  of  these  runs  incorporated  the  Bodner  material 
model  and  the  crack  was  grown  in  the  model  at  a  sufficient  rate  such  that 
model  displacements  matched  test  data  with  time.  Each  node  released  was 
unloaded  in  5  sec.  except  runs  SI  and  S7  as  indicated  in  Table  V-3.  Table 
V-4  summarizes  the  basic  details  for  VISCO  runs  similar  to  Table  V-3  but 
with  no  crack  growth  allowed.  It  can  be  seen  from  these  tables  that  the 
computer  time  required  for  the  high  load  runs  (e.g.  S2  and  A3)  is  much 
higher  than  for  the  lower  K  levels.  In  the  case  of  run  A3,  computer  time 
required  is  high  due  to  a  large  load  causing  a  great  amount  of  plastic 
flow  to  occur.  Recall  from  section  III  that  high  plastic  strain  rates 
result  in  small  time  step  size  which  then  requires  more  times  steps  to  sim¬ 
ulate  a  given  amount  of  time  relative  to  the  case  of  low  plastic  strain 
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Figure  V-lOb. 
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Table  V-3.  Summary  of  Hybrid  Experimental-Numerical  VISCO  Runs 


rates.  Furthermore,  in  the  case  of  run  S2,  extensive  crack  growth  occurred 
requiring  19  nodes  to  be  released.  Each  node  release  also  requires  rela¬ 
tively  small  time  steps  due  to  the  redistribution  of  stress  and  the  associ¬ 
ated  plastic  straining  around  the  cracktip. 

Figures  V-lla  through  V-16a  show  the  match  of  VISCO  displacements  with 
each  particular  test’s  displacement  data.  Also,  the  amount  of  displacement 
in  VISCO  for  no  crack  growth  under  load  is  given.  These  displacements  are 
relative  to  the  displacements  existing  at  the  time  maximum  load  is  achieved 
which  means  all  displacements  prior  to  reaching  maximum  load  were  not 
included  in  this  test  data.  Maximum  load  was  normally  achieved  in  5  sec. 

The  experimental  optical  technique  for  measuring  displacements  is  highly 
sensitive  and  can  resolve  displacements  in  the  neighborhood  of  0.1  micron. 
However  this  technique  loses  sensitivity  when  displacement  measurements  are 
much  larger  than  a  micron.  Since  the  main  interest  was  to  record  displace¬ 
ments  after  reaching  maximum  load  (i.e.  the  creep  crack  growth  data)  and  due 
to  the  desire  to  maintain  high  measurement  resolution  the  load  application 
displacements  were  left  out  of  the  creep  crack  growth  displacement  data. 

Test  6  was  an  exception  where  displacement  measurements  began  after  5  minutes 
of  test  time  had  expired. 

The  main  details  of  the  HEN  procedure  can  be  graphically  seen  in  Fig 
V-lla  and  V-llb.  Observe  in  Fig  V-lla  the  two  specific  curves  of  cracktip 
displacement.  One  from  VISCO  with  no  crack  growth  and  the  other  being  test 
data.  These  curves  start  to  deviate  from  each  other  at  a  time  of  approxi¬ 
mately  4  min.  and  for  a  time  of  20  min.  the  difference  is  fairly  great  as 
is  shown  by  a  bracket  in  Fig  V-lla.  This  bracketed  difference  is  attri¬ 
buted  to  physical  crack  growth  and  requires  release  of  cracktip  nodes  in  this 
simulation. 
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Figure  V-12.  Test  6  HEN  VISCO  Results,  (a)  Displacements  Matched  to  Test  Data 
(b)  Crack  Growth 
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Figure  V-13.  Test  8a  HEN  VISCO  Results,  (a)  Displacements  Matched  to  Test 
Data,  (b)  Crack  Growth 
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Figure  V-16.  Test  12b  HEN  VISCO  Results,  (a)  Displacements 
Matched  to  Test  Data,  (b)  Crack  Growth 


Figures  V-llb  through  V-16b  depict  the  resulting  crack  growth  from 
the  HEN  VISCO  runs  for  each  test.  Note  that  due  to  this  discrete  finite 


element  analysis  technique,  the  crack  growth  is  a  step  funtion  (i.e.  re¬ 
lease  of  individual  cracktip  nodes)  whereas  realistically  the  crack  might 
in  general  be  growing  in  a  smoother  manner  with  time. 

COMPARISON  OF  RESULTS  FROM  USING  DIFFERENT  NODE  UNLOADING  METHODS.  In 

cases  where  the  total  creep  crack  growth  is  only  a  few  node  distances  (7.81  x 
-4 

10  in. -20  microns)  the  displacements  developed  by  VISCO  deviate  signifi¬ 
cantly  from  the  test  displacement  versus  time  curve  (e.g.  see  Fig 
V-15a  which  implies  that  the  model  is  too  compliant  or  the  unloading  of  the 
cracktip  nodes  is  too  rapid.  These  deviations  become  less  significant  for 
larger  amounts  of  crack  growth  as  seen  in  Fig  V-lla.  One  approach  to  mini¬ 
mize  these  deviations  for  small  crack  growth  cases  is  to  make  a  second 
VISCO  run  for  the  same  test.  In  this  second  run,  node  release  times  from 
the  first  VISCO  run  are  input.  Therefore  the  continuous  node  unload  method 
can  then  be  used  as  described  in  section  IV. 

Figure  V-lla  shows  VISCO  displacements  from  employing  the  5  sec.  node 
unload  method  (VISCO  Run  S3)  and  the  continuous  node  unload  method  (VISCO 
Run  S7).  Note  that  the  displacements  from  Run  S7  are  only  slightly  below 
and  run  parallel  to  the  test  data  with  no  excursions.  Also,  Run  SI  in  Fig 
V-13a  matches  test  data  with  no  excursions  and  this  is  a  case  similar  to 
Fig  V-15a  where  overall  crack  growth  was  small  and  deviations  were  large. 

Later  on  crack  growth  criteria  will  be  postulated  based  on  stress  and/ 
or  strain  at  the  cracktip  and  therefore  reliable  data  from  the  HEN  procedure 
is  required  for  these  two  quantities.  Figures  V-17  and  V-18  show  the  differ¬ 
ences  in  stress  and  plastic  strain  at  the  cracktip  which  develop  between 
VISCO  runs  using  the  two  different  node  unload  methods.  For  large  amounts  of 
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Figure  V-17.  Cracktip  Effective  Stress  At  The  Time  The  Cracktip  Node  is  Begun 
To  Be  Released 
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Figure  V-18.  Cracktip  Plastic  Strain  At  The  Time  The  Cracktip  Node  Is  Begun 
To  Be  Released 


crack  growth  in  this  case,  the  relative  difference  in  effective  stress  at 


the  cracktip  is  seen  to  reach  a  steady  state  value  of  approximately  1%.  The 
relative  differences  of  the  y-component  of  plastic  strain  at  the  cracktip, 
as  seen  in  Fig  V-18,  is  much  less.  These  differences  in  stress  and  strain 
are  relatively  small  and  it  is  concluded  that  essentially  the  same  results 
(i.e.  stress  and  strain  at  the  cracktip)  can  be  achieved  for  either  node 
unload  method  for  cases  with  large  amounts  of  crack  growth.  For  crack  growth 
under  small  loads  or  for  small  crack  growth  rates  more  time  is  allowed  for 
stress  relaxation  ahead  of  the  cracktip.  Thus  on  the  average  cracktip 
stress  and  strain  also  differ  little  between  the  two  unload  methods  for  slow 
crack  growth.  With  this  as  background,  it  was  decided  that  all  subsequent 
solutions  would  be  carried  out  using  the  5  sec.  node  unloading  scheme. 

DEPENDENCE  ON  DEFORMATION  HISTORY.  Tests  8b  through  8d  were  not  included 
in  the  HEN  VISCO  runs  since  no  fatigue  precracking  was  done  prior  to  these 
tests.  Without  fatigue  precracking  these  tests  were  considered  to  be  a- 
typical  due  to  their  different  prior  deformation  history.  To  check  out  de¬ 
pendence  on  prior  deformation  history  the  following  VISCO  analysis  was 
performed. 

Three  VISCO  runs  with  no  crack  growth  were  made  to  simulate  Test  8b 
and  its  dependence  on  prior  deformation.  Each  run  had  one  of  three  load¬ 
time  profiles  in  Fig  V-19  input  to  it.  The  VISCO  run  A7  had  no  prior 
load  history,  run  A8  had  Test  8a  load  history  with  a  total  unload  and  re¬ 
load  cycle,  which  is  the  most  realistic  load-time  model  of  Tests  8a  and 
8b,  and  finally  run  A9  had  Test  8a  load  history  with  no  unload,  but  a  load 
increase  to  the  Test  8b  load  level.  The  top  of  Fig  V-19  shows  the  displace¬ 
ment  versus  time  profile  for  each  of  these  VISCO  runs  compared  to  Test  8b 
data.  Again  these  displacements  are  relative  to  the  displacements  existing 
when  the  maximum  Test  8b  load  is  achieved. 
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Figure  V-19.  Load  History  Dependence  of  Displacements 
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In  the  case  of  run  A8  where  complete  unloading  occurs  prior  to  Test 
8b  load,  the  VISCO  results  showed  that  plastic  flow  reoriented  itself  during 
the  unloading  due  to  the  reversing  of  the  principal  stress  at  the  cracktip 
to  compression.  This  stress  develops  during  unloading  since  the  material  has 
prior  tensile  plastic  strains  from  the  Test  8a  load  and  cannot  return  to 
its  original  strain  free  state.  This  compressive  behavior,  if  associated  with 
crack  growth,  leads  to  the  crack  closure  phenomenon  described  by  Newman  [83]. 

Comparing  displacements  of  run  A8  and  A9,  where  no  unloading  was  done, 
indicates  quite  similar  behavior.  When  no  prior  deformation  history  exists 
as  for  run  A7  the  displacements  differ  significantly,  as  shown  in  Fig  V-19 , 
from  those  with  prior  load-deformation  history.  Note  that  run  A7  fits  the 
test  data  curve  best.  However,  test  compliance  data  indicates  that  some 
crack  growth  occurred  in  Test  8b.  If  displacements  associated  with  this 
crack  growth  were  included  in  these  VISCO  runs,  experience  with  Test  8a 
would  indicate  that  runs  A8  and  A9  would  be  brought  up  to  the  test  data  and 
run  A7  would  be  much  in  excess  of  test  results.  It  should  also  be  noted 
that  the  form  of  the  Bodner  model  employed  in  VISCO  does  not  include  the 
Bauschinger  effect  prevalent  in  metals  when  compressive  yielding  follows 
tensile  yielding.  Therefore,  for  the  current  investigation  Tests  8b  through 
8d  will  not  be  further  studied  in  order  to  minimize  dependence  on  pretest  de¬ 
formation  history  in  this  investigation. 

CRACK  GROWTH  RESULTS ■  Table  V-5  displays  the  total  creep  crack  growth 
increments  measured  and  calculated  for  the  indicated  test.  The  column  en¬ 
titled  Aam  presents  the  measured  crack  growth  based  on  the  average  of  four 
measurements  across  the  thickness  as  shown  in  Fig  V-19.  The  column  entitled 
Aa^  is  the  total  crack  growth  calculated  in  the  HEN  VISCO  runs  and  shown  pre- 
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Table  V-5.  Summary  of  Measured  and  Calculated  Crack  Growth 


viously  in  Figs  V-llb  through  V-16b.  Two  other  convenient  methods  for 
calculating  crack  growth  without  the  need  of  HEN  VISCO  runs  were  also  used 


and  their  results  are  presented  as  and  Aa^.  The  method  for  calculating 
Aa^  is  an  original  technique  developed  in  the  present  investigation  and  will 
be  discussed  subsequently. 

The  column  entitled  Aa^  is  the  crack  growth  calculated  by  Clarke's 
elastic  compliance  method  [75].  This  method  predicts  negative  crack  growth 
in  some  instances  where  compliance  was  observed  to  decrease  as  discussed  in 
section  IV.  Figure  V-20  compares  crack  growth  calculated  by  Clarke's  method 
to  the  results  from  a  HEN  VISCO  run  for  Test  9.  Note  that  experimental 
elastic  compliance  changes  from  Test  9  indicate  some  unrealistic  negative 
crack  growth  while  incorporating  Clarke's  method  whereas  the  VISCO  results 
show  a  realistic  monotonically  increasing  amount  of  crack  growth  with  time. 
However,  for  times  greater  than  20  min. ,  both  curves  are  approximately  para¬ 
llel  which  lends  support  to  Clarke's  method  for  large  amounts  of  creep  crack 
growth  also  demonstrated  by  Donat  [76]. 

The  column  in  Table  V-5  entitled  t\a^  is  the  crack  growth  calculated  by 
a  variation  of  Clarke's  elastic  compliance  method  as  described  below.  In 
this  variation  no  unload/reload  cycle  data  are  necessary  during  a  creep  crack 
growth  test  as  discussed  in  section  IV.  The  VISCO  simulation  of  the  test  is 
also  simplified  since  a  VISCO  run  is  made  using  test  conditions,  but  no  crack 
growth  is  allowed  in  the  VISCO  model.  Thus  any  increase  in  displacements 
after  reaching  maximum  load  in  the  no-crack-growth  VISCO  run  can  only  occur 
due  to  time  dependent  plastic  deformation  allowed  by  the  Bodner  material  model. 
Examples  of  these  no-crack-growth  VISCO  displacements  were  given  in  Figs  V-lla 
through  V-16a.  Through  simultaneous  use  of  test  displacement  data  and  the  no- 
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Figure  V-20.  Comparison  of  Crack  Growth  From  Experimental 
Compliance  With  HEN  VISCO  Results 
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crack-growth  VISCO  results,  the  crack  growth  Aa£  is  calculated  as  follows 


C^Tcs4  ~ 

ac/a  a 


(V-6 ) 


where  6  is  the  experimental  displacement.  The  displacement  is  the 
Test  v 

value  calculated  in  VISCO  for  the  same  test  conditions  but  no  crack  growth 
is  allowed.  The  denominator  term,  Ac/Aa  is  the  rate  of  change  of  elastic 
compliance  with  crack  length  generated  from  elastic  VISCO  runs  of  the  test 
specimens.  By  releasing  cracktip  nodes  and  dividing  the  resulting  VISCO 
displacements  by  the  load  to  get  the  respective  compliances,  the  curve  of 
compliance  versus  crack  length  was  determined  for  each  test's  indent  or 
displacement  measurement  location  (see  Fig.  V-6).  The  slope  of  this  com¬ 
pliance  versus  crack  length  curve  then  provided  Ac/Aa  ratio.  The  applied 
load  P  is  divided  into  the  difference  of  6  and  6^  to  generate  a  change  in 
compliance.  This  change  in  compliance  is  assumed  to  be  primarily  dependent 
on  crack  growth  which  means  that  the  plastic  zone  developing  around  this  fixed 
cracktip  is  approximately  the  same  size  for  an  extending  crack  and  simply 
translates  along  with  the  cracktip. 

Figure  V-21  graphically  defines  6  and  6,  at  an  arbitrary  measure- 

Test  V 

ment  time  t  .  Note  that  this  Aan  method  is  effectively  the  same  as  the  HEN 
m  z 

VISCO  procedure  for  small  amounts  of  crack  growth  where,  as  mentioned  above, 
the  plastic  zone  size  of  crack  growth  versus  time  could  also  be  generated  by 
this  Aaj  method  by  applying  Eq  V-6  continuously  along  the  test  data  curve. 

Consider  the  measured  crack  growth  column  Aa^  and  the  HEN  VISCO  results 
Aa^.  The  differences  between  these  two  columns  is  only  +  10%  for  those 
tests  totally  modelled.  Only  part  of  Test  6  was  modelled  since  the  number  of 
uniform  crack  growth  elements  ahead  of  the  initial  cracktip  were  already 
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exhausted  after  half  of  the  test  duration.  This  also  occurred  in  modelling 


Test  9.  Hence,  in  both  Test  6  and  9  the  value  of  Aa^  was  determined  by 
extrapolating  the  final  or  relatively  steady  state  slope  of  the  crack  growth 
versus  time  curves  in  Figs  V-12b  and  V-14b.  Another  special  characteristic 
of  Test  6  is  that  displacement  data  versus  time  was  not  measured  during  the 
first  five  minutes  of  the  test.  Thus,  based  on  other  test  behavior,  a  sig¬ 
nificant  amount  of  the  early  relatively  rapid  displacement  rates  were  ignored. 
Some  of  this  early  displacement  would  apparently  have  required  more  crack 
growth  in  VISCO  than  the  existing  data  did  while  using  the  HEN  procedure. 

The  crack  growth  values,  Aa^,  in  Table  V-5  also  correlate  quite  well 
with  the  measured  values  Aa  .  The  good  correlation  of  Aa„  with  test  data 
provides  support  to  this  convenient  method,  developed  herein,  for  calculating 
creep  crack  growth.  The  significant  convenience  feature  in  calculating  Aa^ 
is  that  once  a  no-crack-growth  VISCO  run  is  made  for  a  given  set  of  test 
conditions,  Eq  V-6  can  simply  be  applied  for  analyzing  all  further  experi¬ 
ments  with  approximately  the  same  test  conditions  (e.g.  load,  crack  length 
and  geometry) . 

A  somewhat  similar  HEN  analysis  was  done  by  HSU  et.  al.  [32] to  predict 
crack  growth  in  zirconium  at  elevated  temperatures.  HSU's  crack  growth  pre¬ 
dictions  were  2.5  times  greater  than  actual  test  data.  In  the  present  HEN 
analysis,  crack  growth  predictions  were  within  approximately  10%  of  test  data 
for  the  tests  that  were  totally  simulated.  The  extremely  good  correlation  using 
the  present  method  is  attributed  to  both  higher  resolution  experimental  dis¬ 
placement  data  than  HSU's  and  a  more  realistic  material  model  that  includes 
creep  behavior.  Additional  displacements  due  to  creep  in  HSU's  analysis  would 
have  reduced  the  amount  of  crack  growth  predicted.  Since  HSU's  predictions  were 
high,  this  reduction  would  be  in  the  direction  of  better  correlation  with  test 
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Another  comparison  of  current  results  can  qualitatively  be  made  with 
Newman’s  finite  element  analysis  of  fatigue  crack  propagation  [83].  Figure 
V-22  shows  crack  opening  displacement  profiles  during  creep  crack  growth 
from  the  VISCO  model  as  represented  by  the  solid  lines.  Note  that  with  the 
cracktip  at  point  zero,  the  elastic  crack  opening  displacement  profile (i.e. 
the  dashed  line  in  Fig  V-22)  is  lower  than  the  solid  line  which  incorporates 
the  Bodner  material  model.  However,  as  the  crack  grows  a  wake  of  residual 
plastic  deformation  is  left  behind  the  cracktip.  After  12  increments  of 
crack  growth.  Fig  V-22  shows  how  this  plastic  wake  diminishes  the  elastic- 
plastic  crack  profile  below  a  purely  elastic  profile  for  the  same  crack  length 
"a".  The  residual  plastic  deformation  indicated  by  the  cross-hatched  area 
in  Fig  V-22  was  also  displayed  in  a  similar  figure  by  Newman. 

CRACK  GROWTH  CRITERION 

Based  on  the  good  correlation  between  actual  and  predicted  crack  growth, 
the  VISCO  results  from  the  HEN  applications  were  examined  for  potential  crack 
growth  parameters.  The  local  cracktip  parameters  such  as  strain  and  crack 
opening  displacement  (C.O.D.)  were  examined  initially  since  these  parameters 
have  shown  promise  elsewhere  [84,  85]  for  correlating  finite  element  results 
with  crack  growth  test  data. 

The  main  goals  here  were  to  check  out  the  validity  of  existing  crack 
growth  criteria  (i.e.  critical  C.O.D.  and  strain),  possibly  modify  one  of 
these  existing  criteria  to  better  fit  test  data,  and  postulate  a  new  criterion 
that  might  better  account  for  creep  damage  accumulation  and  crack  growth  dis¬ 
played  by  the  HEN  VISCO  analysis. 

CRITICAL  STRAIN  CRITERION.  Examination  of  the  strains  in  the  elements 
adjacent  to  the  cracktip  prior  to  each  node  release  for  crack  growth  in  the 
HEN  VISCO  runs  revealed  that  no  single  value  for  the  critical  strain  would 
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Elastic  VISCO  Results 


Figure  V-22.  Crack  Opening  Displacement  Profiles 


satisfy  all  test  conditions.  However,  a  few  VISCO  runs  were  made  using 
a  fixed  critical  strain  criterion  to  further  evaluate  its  applicability. 

The  critical  strain  criterion  was  implemented  in  VISCO  by  comparing 
the  average  of  the  plastic  strain  components  normal  to  the  crack  and  with¬ 
in  elements  adjacent  to  the  cracktip  with  a  critical  strain  value,  ecrit> 
as  time  progressed.  When  the  average  cracktip  plastic  strain  exceeded  the 
critical  value,  Ecr^t>  the  cracktip  node  was  released  and  unloaded  in  5  sec. 
Table  V-6  gives  a  summary  of  the  basic  details  for  VISCO  runs  employing  the 

critical  strain  criterion.  A  e  .  value  of  0.030  was  found  to  work  well 

crit 

from  the  HEN  VISCO  runs  of  Test  8a  where  was  16.3  ksi/in.  However  for 

Test  9  where  was  36.8  ksiyiii  the  ecrit  value  needed  to  be  0.090  to  work 
well  as  shown  in  Fig  V-23a.  The  corresponding  ecrit  dependent  crack  growth 
versus  time  is  given  in  Fig  V-23b.  The  ecr^t  value  of  0.075  allowed  too 
much  displacement  compared  to  test  data  as  seen  in  Fig  V-23.  The  displace¬ 
ment  results  for  a  value  of  0.090  appear  to  fit  the  data  quite  well  over 

the  first  five  min.  of  the  test.  Likewise  the  resulting  crack  growth  in  Fig 
V-23b  agrees  quite  well  with  the  HEN  data  in  Fig  V-13b.  However,  based  on 
examination  of  the  prior  HEN  VISCO  runs,  if  the  simulation  time  were  continued 
the  Ecrit  value  of  0.090  would  have  been  too  large.  Therefore,  insufficient 
crack  growth  would  have  resulted  and  the  displacements  would  have  fallen  away 
from  the  test  data  as  shown  by  the  extrapolated  dashed  curve. 

Figure  V-24  displays  the  HEN  VISCO  cracktip  strain  values  taken  at  the 
time  a  cracktip  node  was  to  be  released  in  the  HEN  VISCO  runs.  This  plot  was 
motivated  in  an  attempt  to  find  a  general  trend  of  the  critical  strain  values 
or  to  determine  a  mathematical  relationship  with  time  to  envelope  the  HEN 
results.  Note  that  in  general  the  strains  are  high  for  short  times  and  then 
diminish  with  time  to  a  common  value  of  approximately  0.03.  In  order  to  de- 
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Figure  V-2A.  Time  Dependent  Critical  Strain  Crack  Growth  Criterion 
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velop  a  critical  strain  functional  expression  one  can  see  a  need  for 
a  decaying  parameter  which  would  fit  the  upper  bound  strain  values  in 

Fig  V-24  related  to  rapid  crack  growth.  In  addition  the  expression 
must  have  the  capability  of  allowing  the  region  around  the  cracktip 
node  to  deteriorate  with  exposure  to  the  cracktip  environment  (i.e.  not 
only  the  environment  external  to  the  specimen  but  also  the  singular  strain 
state  internal  to  the  specimen) .  Cracktip  deterioration  is  considered  to 
be  displayed  by  the  lower  bound  of  strain  values  in  Fig  V-24.  Therefore 
the  general  properties  of  the  critical  strain  function  have  been  stated. 
The  curve  must  be  decaying  with  test  time,  thus  a  negative  exponential 
function  is  in  order.  And  since  critical  strain  values  appear  to  diminish 
with  cracktip  exposure  time,  but  not  as  rapidly  as  an  exponential  function 
would  dictate,  the  cosine  function  was  examined. 

An  empirical  expression  for  the  critical  strain  that  fits  the  HEN 
VISCO  results  fairly  well  is  given  by 


£cr'.+ 


Co  [  A  exp(-bt)  COS  (££)  +  j ) 
6* 


if  T  -  T0 

if  T^re 


( V  —  7 ) 


where:  e  =0.03,  A=3.0,  b=l . 34xl0_3sec_1 ,  T  =600. sec. 
o  o 

The  value  of  eq  represents  the  critical  strain  for  large  time.  The  co¬ 
efficient  A  is  determined  at  t  equal  to  zero.  Once  A  and  are  chosen 
b  is  determined  by  best  fitting  the  upper  bound  where  cracktip  exposure 

time  is  small  and  set  to  zero  (i.e.  T=o) .  The  parameter  T  determines  how 

o 

rapidly  the  critical  strain  diminishes  to  gq  with  cracktip  exposure  time 


T. 


Motivation  for  the  development  of  Eq  V-7  comes  from  environmental 
effects  such  as  oxidation  at  these  high  temperatures  from  exposure  to 
laboratory  air.  This  oxidation  is  then  associated  with  changing  material 
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properties  at  the  cracktip  such  as  the  critical  strain  value.  Crack  growth 
in  alloys  similar  to  IN-100  has  been  found  to  be  quite  sensitive  to  the  en¬ 
vironment  and  the  resulting  oxidation  that  can  occur  when  an  elevated  tem¬ 
perature  crack  growth  test  is  done  in  air  [86].  Therefore  Eq.  V-7  is  an 
attempt  to  represent  the  rate  at  which  the  critical  strain  for  crack  growth 
is  diminished  with  time  due  to  environmental  effects. 

The  critical  strain  criterion  just  formulated  would  be  employed  in 
VISCO  in  the  following  steps. 

1.  Register  both  total  test  simulation  time,  t,  and  cracktip  exposure 
time,  T. 

2.  Evaluate  the  critical  strain  level  during  each  time  step  from  Eq . 
V-7. 

3.  Compare  with  cracktip  element’s  plastic  strain  accumulating  in 
the  VISCO  analysis. 

4.  Release  cracktip  node  when  VISCO  plastic  strain  at  the  cracktip 
exceeds  Ecrit  from  Eq .  V-7. 

5.  Cracktip  exposure  time  T  is  set  to  zero  and  the  above  steps  1-4 
are  repeated  for  the  next  cracktip  node. 

CRITICAL  CRACK  OPENING  DISPLACEMENT  CRITERION.  The  crack  opening  dis¬ 
placement  is  here  defined  as  the  C.O.D.  at  the  first  node  behind  the  crack¬ 
tip  in  the  finite  element  model.  Examination  of  C.O.D. 's  taken  from  HEN 
VISCO  runs  in  Table  V-3  just  prior  to  releasing  a  cracktip  node  revealed 
no  single  C.O.D.  value  that  could  be  used  for  all  test  conditions  as  a  cri¬ 
tical  C.O.D.  A  value  of  0.280  x  10  ^  in.  was  found  from  HEN  VISCO  results 
for  Test  8a  with  K^j.  "16.3  ksi^Tn  to  work  best,  yet  when  a  few  VISCO  runs, 
as  summarized  in  Table  V-7,  were  done  with  a  critical  C.O.D.  criterion  and 
an  increase  to  1^=36. 8  ksi)H!n  for  Test  9,  the  best  critical  C.O.D.  became 
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It  is  recognized  that  Eq  V-9  applies  most  directly  to  creep  rupture 
of  uncracked  uniaxial  components.  However,  the  present  author  considers 


that  similar  behavior  might  be  possible  in  creep  crack  propagation.  A 
schematic  of  the  postulated  behavior  involved  with  creep  crack  propagation 
is  given  in  Fig  V-26 .  This  figure  shows  a  creep  damage  front  preceding  the 
crack.  Within  this  front,  the  material  is  accumulating  creep  damage  in  the 
form  of  microcracks.  This  type  of  creep  damage  is  also  associated  with 
creep  rupture  of  uncracked  components. 

In  order  to  apply  Eq  V-10  as  a  crack  growth  criterion  a  process  zone, 

<5,  is  required  and  defined  in  Fig  V-26.  In  addition,  the  rupture  time 
t  is  redefined  from  Eq  V-8  to  the  elapsed  time  the  crack  requires  to  grow 
from  one  node,  in  the  HEN  VISCO  results,  to  the  next  node.  In  other  words, 
it  is  the  time  period  during  which  the  process  zone,  6  is  exposed  to  the 
cracktip  stress  field  prior  to  rupture.  In  the  VISCO  finite  element  analysis, 
this  process  zone  was  taken  as  one  element  preceding  the  cracktip.  The 
average  component  of  stress  normal  to  the  crackline  from  three  elements 
adjacent  to  the  cracktip  was  used  as  the  stress,  6,  in  Eq  V-10  as  shown  in 
Fig  V-26. 

Since  the  greatest  stress  exists  at  the  cracktip  and  environmental  de¬ 


gradation  is  considered  to  be  most  prevalent  there,  most  damage  accumulation 
was  assumed  to  occur  in  the  process  zone  after  the  arrival  of  the  cracktip 
to  the  process  zone's  border.  Therefore  time  in  Eq  V-10  was  measured  from 
cracktip  arrival  time  to  the  current  cracktip  node,  t  ,  or 

,Vtr 

/  o~M  d  t  =  C 

rv-i 


(V-ll) 


The  constants  M  and  C  were  determined  based  on  results  from  the  HEN  VISCO 


runs.  To  accomplish  this,  Eq  V-ll  was  approximated  as 


( V-12) 


The  rupture  time  or  crack  growth  times  t  were  taken  from  HEN  V1SC0  results 
The  stress,  °aVg’  also  based  on  HEN  results,  was  an  average  over  time  t 
of  the  cracktip  stress  defined  previously.  Since  the  interest  here  is  to 
develop  values  for  M  and  C  which  apply  to  the  entire  set  of  tests,  it  be¬ 
came  obvious  that  there  were  more  combinations  of  t  and  oavg  than  nece¬ 
ssary  to  uniquely  define  ..  .  nd  C.  Consequently,  to  include  the  data  from 
each  HEN  VISCO  run,  a  le.  .  ->.  uare  fit  of  the  data  on  a  log-log  plot  was 
used  to  best  fit  the  time-stress  data.  The  values  chosen  were 


M=15 

C-8.63  x  1079  (psi)^5  sec 

Equation  V-ll  along  with  the  above  constants  were  then  incorporated 
into  VISCO  as  a  critical  damage  accumulation  criterion  for  crack  growth. 

Table  V-8  summarizes  the  basic  details  of  the  VISCO  runs  using  this  crack 
growth  criterion. 

Figures  V-27  through  V-32  give  the  results  of  applying  VISCO  with  the 
critical  damage  accumulation  criterion  to  the  indicated  test  number  in  each 
figure.  The  resulting  VISCO  displacements  are  compared  to  test  data  in  Figs 
V-27a  through  V-32a  whereas  the  resulting  crack  growth  from  VISCO  is  given 
in  Figs  V-27b  through  V-32b.  For  the  higher  loadings  such  as  Test  9  in  Fig 
V-30a  the  difference  between  test  data  and  VISCO  results  is  greatest.  Part 
of  this  difference  may  be  due  to  the  fact  that  at  this  high  load  (Kj.£=36.8 
ksiyTn)  the  damage  zone  as  indicated  in  Fig  V-26  is  much  larger  than  the  pro¬ 
cess  zone  sited  used  herein  (i.e.  one  element  size  or  a  characteristic  dimen- 
-4 

sion  of  7.81  x  10  in.).  Hence  significant  damage  accumulation  may  occur 

in  the  material  before  arrival  of  the  cracktip  or  time  t^  for  a  given  element 
in  the  crackpath.  This  damage  occurring  in  an  element  prior  to  t  for  that 
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CRACK  GROTH  (nodes)  TIP  DISPLACEMENT  (microns) 


TIME  (min) 


Figure  V-27.  Test  5  -  VISCO  Critical  Damage  Accumulation  Criterion, 
(a)  Displacement  Results,  (b)  Crack.  Growth  Results 


CRACK  GROWTH  (nodes)  TIP  DlSPLACEMENT(microns) 


Figure  V-28.  Test  6  -  VISCO  Critical  Damage  Accumulation  Criterion, 
(a)  Displacement  Results,  (b)  Crack  Growth  Results 
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CRACK  GROWTH  (nodes)  TIP  DISPLACEMENT  (microns) 


TIME  (min) 


Figure  V-29.  Test  8a  -  VISCO  Critical  Damage  Accumulation  Criterion, 
(a)  Displacement  Results,  (b)  Crack  Growth  Results 
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Figure  V-32. 


Test  12b  -  VISCO  Critical  Damage  Accumulation 
Criterion,  (a)  Displacement  Results,  (b)  Crack 
Growth  Results 
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respective  element  was  neglected  in  the  present  damage  accumulation  criterion. 
However,  for  the  lower  loads  agreement  with  the  test  data  is  quite  good 
considering  that  creep  rates  for  the  same  test  conditions  can  easily 
vary  by  a  factor  of  two  or  three  which  is  why  log-log  plots  are  used  to 
plot  creep  data  [29]. 

It  should  also  be  noted  that  the  M  and  C  values  used  were  determined 
from  the  approximate  Eq  V-12  expression.  Crack,  growth  results  from  this 
criterion  might  also  be  improved  by  iterating  or  making  small  changes 
to  M  and  C  and  making  further  VISCO  runs  in  an  effort  to  better  fit  test 
data. 


CRACK  GROWTH  RATE  CRITERIA 

In  this  section  crack  growth  rate  criteria  will  be  discussed  based  on 
the  steady  state  crack  growth  rates  developed  by  the  HEN  VISCO  runs.  The 
steady  state  crack  growth  rates  were  determined  by  best  fitting  a  straight 
line  to  the  crack  growth  curves  in  Figs  V-llb  through  V-16b.  The  initial 
transient  portion  of  these  crack  growth  curves  was  ignored  for  these  steady 
state  values.  Also  a  discussion  of  the  merits  of  the  C*  integral  as  a 
crack  growth  rate  criterion  will  be  given. 

STRESS  INTENSITY  FACTOR  CRITERION  .  Crack  growth  rate,  a,  has 
previously  been  found  to  correlate  with  the  elastic  stress  intensity 
factor  as  given  in  Eq  II-2  and  recalled  here 


This  relationship  plots  as  a  straight  line  on  log-log  paper,  as  shown  in 
Fig  V-33.  The  experimental  data  referred  to  in  Fig  V-33  was  for  IN-100 
behavior  at  1350°F  which  is  the  same  alloy  and  temperature  used  in  the 
present  investigation.  Note  that  Donat's  experimental  data  [76]  covered 
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CRACK  GROWTH  RATE  (in/hr) 


STRESS  INTENSITY  FACTOR  K  (ksi  i^n) 

Figure  V-33.  Crack  Growth  Rate  vs.  Stress  Intensity  Factor 
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a  range  in  K  values  from  30  to  approximately  80  ksi  in.  In  order  to 
compare  with  the  lower  K  levels  in  the  current  investigation,  the  line 
representing  the  best  fit  to  Donat's  data  was  extrapolated  as  shown  by 
the  dashed  line  in  Fig  V-33.  Agreement  with  the  present  HEN  VISCO  re¬ 
sults,  in  which  Kpg  is  taken  from  Table  V-3  and  a  from  Figs  V~llb  through 
V-16b  is  good  especially  considering  the  fact  that  the  test  data  line 
was  extrapolated. 

This  criterion  has  the  distinct  advantage  relative  to  criteria  pre¬ 
sented  earlier  that  once  the  constants  A  and  a  are  determined,  it  can  be 
used  independent  of  finite  element  analyses.  This  advantage  is  due  to  the 
fact  that  K  can  be  calculated  for  most  test  geometries  by  relatively  simple 
equations  like  Eq  V-3.  Thus  the  so-called  steady  state  creep  crack  growth 
rate  can  be  simply  calculated  from  Eq  V-13,  but  it  should  be  kept  in  mind 
that  incubation  time  for  crack  initiation  nor  the  initial  rapid  crack  growth 
observed  in  the  HEN  VISCO  results  is  captured  with  this  criterion. 

NET  SECTION  STRESS  CRITERION.  Crack  growth  rates  have  also  been  shown 


to  be  related  to  net  section  stress  as  given  in  Eq  II-3  and  recalled  here 

a  =■  B  (cr^  )  ( V— 14) 


Figure  V-34  shows  how  Eq  V-13  also  plots  as  a  straight  line  (the  dashed  line) 
on  log-log  paper.  The  present  results  seem  to  correlate  to  Eq  V-14  as  well 
as  they  did  to  Eq  V-12  on  log-log  paper.  This  is  not  surprising  since  for 
the  center  cracked  plate  it  can  be  shown  that  K  is  approximately  directly 
proportional  to  for  crack  lengths  of  aA\’  from  .2  to  .7  which  spans  the 
the  crack  lengths  in  the  current  study.  The  net  section  stress  criterion, 
like  the  K  criterion,  also  neglects  crack  growth  characteristics . 

C*  INTEGRAL  CRITERION.  The  C*  integral  is  an  extension  of  the  J  integral 


concept  for  application  to  creep  crack  growth  (a  detailed  discussion  of  C* 
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CRACK  GROWTH  RATE  a  (in/hr) 


Figure  V-34.  Crack  Growth  Rate  vs.  Net  Section  Stress 
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is  given  in  Appendix  C) .  Theoretically  the  C*  integral  is  path  independent 
for  a  creeping  solid  where  stress  is  only  a  function  of  the  plastic  strain 
rate  and  elastic  strain  rates  do  not  enter  into  the  formulation.  In  the 
present  research,  C*  was  calculated  as  shown  in  Appendix  C.  In  this  cal¬ 
culation  total  strain  rates  were  attributed  to  creeping  plastic  strain 
rates  which  implies  no  elastic  strain  rates  exist.  For  a  realistic  materia] 
that  includes  elastic  behavior,  if  the  elastic  strain  rates  are  zero  then 
the  stress  state  is  constant  with  time.  And  therefore,  the  creeping  plas¬ 
tic  strains,  since  they  are  a  function  of  stress  must  also  be  constant. 
Accordingly,  with  constant  stress  and  strain  rate  the  W*  integral  was  in¬ 
tegrated  directly  in  Appendix  C. 

The  above  description  for  calculating  C*  in  VISCO  was  implemented 
and  applied  to  a  realistic  elastic-viscoplastic  material  (i.e.  IN-100) 
being  studied  herein.  It  was  observed  that  during  load  application 
the  C*  values  were  extremely  high  due  to  the  elastic  strain  rate  con¬ 
tribution.  After  maximum  load  was  achieved  C*  values  reduced  down  to 
much  lower  steady  state  values  until  crack  growth  began  in  the  HEN  VISCO 
runs.  Again  as  crack  growth  began,  the  C*  values  increased  significantly 
due  to  the  elastic  strain  rate  contribution  as  stress  was  redistributing 
around  the  moving  cracktip.  It  should  be  noted  that  any  attempt  to  re¬ 
move  elastic  contributions  to  strain  and  displacements  will  result  in  an 
ill  posed  problem  since  the  displacement  rate  components  needed  in  Eq  11-7 
cannot  be  resolved  into  an  elastic  and  plastic  portion. 

Thus  far  evaluation  of  the  C*  integral  given  in  Eq  1 1  —  7 ,  the  elastic 
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part  of  the  strain  rate  is  neglected,  otherwise  the  integral  has  no 
meaning.  But  when  the  elastic  strain  rates  are  ignored,  C*  should  only 
be  calculated  after  the  stresses  are  fairly  constant  and  prior  to  crack 
initiation,  C*  will  then  be  a  constant  until  the  crack  begins  to  grow. 
Therefore  it  is  impossible  to  relate  C*  to  any  incubation  time  for  crack 
growth.  Moreover,  during  crack  growth  the  contribution  from  elastic  strain 
rates  is  again  substantial  and  cannot  be  ignored.  Therefore  it  appears  that 
the  C*  integral  is  ineffective  as  a  fracture  criterion  in  a  finite  element 
model  for  creep  crack  growth  in  the  current  investigation. 

LOAD  POINT  DISPLACEMENT  RATE  CRITERION.  An  expression  relating  crack 
growth  rates  to  load  point  displacement  rate  was  given  in  Eq  II-4.  Unfor¬ 
tunately  this  criterion  suffers  from  problems  similar  to  the  C*  integral. 

The  load  point  displacement  rate  after  reaching  maximum  load  is  the  sum  of 
the  displacement  rate  due  to  crack  growth  as  well  as  the  displacement  rate 
due  to  plastic  deformation.  Hence  the  load  point  displacement  rate  may  have 
several  values  for  the  same  crack  growth  rate  depending  on  the  rate  of 
plastic  deformation  such  as  demonstrated  in  Eq  V-6.  One  can  get  an  appre¬ 
ciation  of  how  much  plastic  deformation  contributes  to  the  overall  displace¬ 
ments  by  looking  at  the  no-crack-growth  displacements  versus  time  in  Figs 
V-lla  through  V-16a  which  are  totally  a  result  of  plastic  deformation.  There¬ 
fore  the  load  displacement  rate  does  not  provide  a  unique  solution  to  the 
crack  growth  rate  unless  variations  in  the  plastic  deformation  rate  can  be 
neglected. 
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VI  SUMMARY  AND  CONCLUSIONS 


A  two  dimensional  (plane  stress/plane  strain)  finite  element  program 
has  been  developed  which  accounts  for  both  nonlinear  viscoplastic  material 
behavior  and  changing  boundary  conditions  due  to  crack  growth.  Three  visco¬ 
plastic  material  models:  (1)  Malvern  Flow  Law,  (2)  Norton’s  Creep  Law,  and 
(3)  Bodner-Partom  Flow  Law  were  incorporated  into  the  program.  These  time 
dependent  material  models  were  numerically  integrated  through  time  by  a 
linear  Euler  extrapolation  technique.  A  variable  time  step  algorithm  was 
included  that  maximized  time  step  size  during  the  analysis  while  maintain¬ 
ing  good  accuracy.  This  program  was  used  as  the  plane  stress  theoretical 
model  for  the  hybrid  experimental-numerical  procedure  employed  to  analyze 
sustained  load  creep  crack  growth  test  data.  The  test  specimens  were 
center  cracked  plates  made  of  IN-100  and  tested  at  1350°F. 

The  following  statements  and  conclusions  are  based  on  the  creep  crack 
growth  analysis  herein. 

1.  A  method  for  getting  crack  growth  behavior  solely  from  displacement 
measurements  in  conjunction  with  a  cracked  specimen  model  which  uti¬ 
lizes  realistic  constitutive  relationships  has  been  developed.  The 
constitutive  law  was  especially  tailored  to  the  nickel-base  alloy 
studied  which  displays  time  dependent  nonlinear  inelastic  behavior 
at  elevated  temperatures.  It  has  been  demonstrated  that  the  technique 
can  be  applied  where  crack  extension  is  very  small  and  could  not 
otherwise  be  resolved  by  conventional  experimental  crack  measuring 
techniques  (e.g.  compliance  techniques  or  using  a  travelling  micro¬ 
scope).  This  method  provides  realistic  monotonically  increasing 
crack  growth  values  with  resolution  better  than  0.001  in.  Crack  growth 
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2. 


predictions  agreed  to  within  10%  of  post-test  measurements. 

Numerical  procedures  were  developed  to  efficiently  integrate  the  non¬ 
linear  time  dependent  material  models  and  simulate  crack  growth  by  two 
cracktip  node  unloading  methods.  The  numerical  time  integration  pro¬ 
cedure  utilized  an  Euler  linear  extrapolation  technique  with  a  vari¬ 
able  time  step  algorithm  that  maximized  time  step  size  during  the 
simulation  while  maintaining  good  accuracy.  One  cracktip  node  unload 
method  diminished  the  restraining  force  on  the  finite  element  cracktip 
node  in  a  specified  5  sec.  time  period  independent  of  crack  growth  rate 
whereas  the  second  method  continuously  unloaded  cracktip  nodes  in  pro¬ 
portion  to  a  predetermined  crack  growth  rate.  Unloading  the  nodes  con¬ 
tinuously  provided  a  closer  fit  to  displacement  versus  time  test  data 
however  the  average  displacement  versus  time  was  approximately  the  same 
for  both  node  unloading  methods. 

3.  A  procedure  was  developed  for  determining  crack  extension  using  cal¬ 
culations  of  viscoplastic  deformation  with  no  crack  growth.  In  this 
procedure,  the  difference  between  total  test  deformation  and  viscoplastic 
deformation  is  attributed  to  crack  extension.  Extremely  good  crack 
growth  predictions  were  made. 

A.  The  elastic  compliance  method  for  resolving  creep  crack  extension  has 

been  shown  to  imply  negative  unrealistic  crack  growth  and  is  unreliable 
especially  during  the  first  part  of  a  creep  crack  growth  test. 

5.  Several  parameters  were  studied  for  their  potential  as  creep  crack  growth 
controlling  parameters. 

a.  No  single  fixed  value  of  strain  for  a  critical  strain  crack  growth 

criterion  was  found  to  match  all  test  conditions  in  this  investigation. 
Environmental  effects  apparently  tend  to  lower  the  critical  strain 
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magnitude  with  time,  under  load.  An  empirxcal  relationship  was 
developed,  based  on  the  HEN  results,  which  gives  the  critical 
strain  a  diminishing  value  with  time.  This  investigation  did  not 
include  any  plane  strain  analyses.  Due  to  higher  constraint  at 
the  cracktip  for  plane  strain,  less  stress  redistribution  would 
occur.  Therefore  it  seems  possible  that  the  critical  strain  for 
all  test  conditions  may  vary  less  in  a  plane  strain  simulation. 

b.  No  single  fixed  value  for  C.O.D.  was  found  to  match  all  test  con¬ 
ditions  using  a  critical  C.O.D.  crack  growth  criterion.  The  C.O.D. 
behaved  similar  to  the  cracktip  strain  with  time  however  its  per¬ 
cent  variation  was  less. 

c.  A  critical  damage  accumulation  criterion  for  crack  growth  was  de¬ 
veloped  based  on  a  modification  of  the  life  fraction  rule  for  creep 
rupture  to  account  for  environmental  effects  at  the  cracktip.  App¬ 
lication  of  this  criterion  provided  good  agreement  with  the  low  to 
medium  load  test  conditions.  For  the  highest  load  test  cases,  this 
criterion  predicted  crack  growth  rates  somewhat  lower  than  the  HEN 
results.  It  appears  that  accumulation  of  damage  over  all  time  and 
not  just  cracktip  exposure  time  might  improve  the  results. 

Data  obtained  in  this  Investigation  through  numerical  calculations  pro¬ 
vided  crack  growth  versus  time  and  thus  crack  growth  rate  time,  a  . 

The  a  data  compared  well  with  published  data  for  the  same  material  and 
temperature  when  plotted  against  stress  intensity  factor.  The  present 
data  was  obtained  for  a  and  K  values  lower  than  the  referenced  data. 

Net  section  stress  also  provided  good  correlation  with  the  predicted 
crack  growth  rates. 

The  C*  integral  and  load  line  displacement  rate  were  investigated  as 


possible  parameters  controlling  crack  growth  rate,  a.  The  C*  in¬ 
tegral  is  an  unreliable  parameter  for  predicting  creep  crack  growth 
due  to  its  formulation  which  is  based  on  a  creeping  solid  behavior 
that  neglects  elastic  strain  rates.  The  load  line  displacement  rate 
which  can  be  shown  to  be  proportional  to  C*  also  does  not  provide  a 
unique  solution  for  the  crack  growth  rate  unless  variations  in  plastic 
deformation  rate  can  be  ignored.  In  general  these  parameters  appear  to 
have  no  applicability  to  crack  growth  rate  prediction  using  numerical 
modeling  of  materials.  However  these  parameters  seem  to  correlate  a 
data  fairly  well  once  the  solution  is  known  as  seen  in  the  literature. 

8.  Constant  strain  triangular  finite  elements  of  the  size  of  grains  at  the 
cracktip  work  well  for  resolving  small  increments  of  creep  crack  growth 
through  the  method  developed  herein. 

9.  VISCO  results  using  the  Bodner-Partom  material  model  were  very  similar 

to  the  VISCO  results  incorporating  the  Malvem-Norton  superposition  model 
for  the  center  cracked  plate,  especially  for  times  greater  than  200  sec. 
after  load  application.  Also  for  times  greater  than  200  sec.,  Norton's 
law  alone  in  VISCO  was  very  similar  to  VISCO  results  using  the  Bodner- 
Parton  material  model. 

10.  Displacements  and  the  associated  crack  growth  were  found  to  be  signifi¬ 
cantly  dependent  on  prior  deformation  history.  Prior  deformation  history 
became  very  important  in  the  case  of  Test  8  where  no  fatigue  precracking 
was  done  between  creep  crack  growth  tests. 

The  above  advancements  in  the  understanding  of  creep  crack  growth  behavior 
at  elevated  temperature  are  especially  suited  for  aiding  future  slow  crack  growth 
tests  for  determining  the  threshold  load  levels  for  creep  crack  growth.  In 
addition  the  crack  growth  criteria  investigations  provide  significant  progress 
towards  life  prediction  of  actual  turbine  disks. 
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The  following  additional  research  work  is  recommened  to  further  the 


life  prediction  capability  developed  in  the  present  work. 

1.  Use  present  approach  to  analyze  other  test  specimen  geometries  (e.g. 
compact  tension  specimen)  to  determine  dependence/ independence  of  results  on 
specimen  geometry  and  also  assess  repeatability  and  material  data  scatter. 

2.  As  computers  become  faster  and  more  efficient  finite  element 
techniques  are  developed,  a  three  dimensional  analysis  of  creep  crack  growth 
should  be  accomplished  to  correctly  model  through  the  thickness  variat ions 
in  test  specimen  behavior.  As  a  first  step  in  this  direction,  plane  strain 
analyses  similar  to  the  present  work  might  provide  additional  insight  into 
creep  crack  growth  behavior. 

3.  Future  work  needs  to  include  cyclic  or  engine  spectrum  loading  conditions. 

4.  Additional  material  characterization  test  data  is  needed  in  general  for 
IN-100.  Tests  providing  this  data  should  be  done  at  several  temperatures  such 
that  the  constitutive  model  could  be  further  developed  and  include  temperature 
dependence . 

5.  Environmental  effects  should  be  further  researched  by  performing  creep 
crack  growth  testing  in  several  different  environments  such  as  vacuum,  inert, 
salt  spray,  high  sulfur  content,  variable  oxygen  partial  pressures,  etc. 

In  summary  the  technique  developed  herein  is  worth  exploring  further  in  that 
it  has  potential  for  providing  information  on  crack  growth  rate  behavior  in 
engine  materials  under  typical  operating  conditions  which  might  not  be  readily 
obtained  using  conventional  techniques.  This  information  in  turn,  is  necessary 
in  order  to  implement  a  retirement  for  cause  philosophy  for  U.S.  Air  Force  jet 
engine  components. 
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APPENDIX  A.  FINITE  ELEMENT  FORMULATION 


In  the  following  sections,  the  basic  concepts  and  equations  used  in 
the  finite-element  analysis  of  elastic  and  elastic-plastic,  materials  are 
briefly  reviewed.  These  equations  provide  background  for  the  elastic- 
viscoplastic  nonlinear  finite-element  computer  program  development. 

The  basic  philosophy  of  the  finite-element  method  [31]  is  that  an 
approximate  solution  to  a  complicated  problem  can  be  obtained  by  sub¬ 
dividing  the  region  of  interest  into  a  finite  number  of  discrete  elements 
and  then  choosing  appropriate  relatively  simple  functions  to  represent 
the  solution  within  each  element.  These  functions  are  simple  compared 
to  the  so-called  "exact"  solutions  which  account  for  the  entire  region  of 
interest.  In  this  section  the  equations  associated  with  representing  a  two- 
dimensional  body  as  a  finite  number  of  elements  are  presented.  The  dis¬ 
placements  in  each  element  were  expressed  as  a  simple  polynomial  and  the 
equations  relating  displacements  to  applied  loading  for  both  plane-stress 
and  plane-strain  conditions  are  given. 

DISPLACEMENT  MODEL 

The  displacement  function  used  in  the  displacement  formulation  is 
generally  selected  as  a  polynomial.  The  polynomial  expression  allows  for 
simple  differentiation  and  integration.  Also,  as  the  element  size  becomes 
small,  the  polynomial  expression  permits  a  simple  approximation  to  the 
exact  solution.  A  polynomial  of  infinite  order  corresponds  to  an  exact 
solution.  However,  for  practical  purposes  the  polynomial  must  be  truncated 
to  a  finite  number  of  terms.  Thus,  the  number  of  elements  in  a  structure 
must  be  large  enough  so  that  the  displacement  function  for  each  element 
closely  approximates  the  exact  displacements  in  that  particular  region. 

In  any  numerical  method,  the  solution  should  converge  to  the  exact 


150 


solution  as  the  size  of  the  elements  become  small.  For  the  displacement 
formulation,  it  has  been  shown  that  under  certain  conditions  Lite  solution 
provides  a  lower  bound  to  the  exact  displacements  [31].  To  assure  this 
convergence  certain  conditions  must  be  satisfied.  First,  the  displacement 
function  must  be  chosen  so  that  rigid  body  displacements  do  not  cause 
straining  of  the  element.  Second,  the  function  must  also  be  chosen  so  that 
a  constant  state  of  strain  is  obtained  as  the  element  size  approaches  zero. 

The  simplest  polynomial  function  which  staisfies  these  two  requirements  and 
also  maintains  displacement  continuity  between  adjacent  elements  is  the 
linear-displacement  function. 

Displacement  Function.  Figure  A-l  shows  a  typical  triangular  element,  m, 
with  nodes  i,  j,  k  numbered  in  a  counter-clockwise  direction.  The  linear- 
displacement  function  which  defines  the  displacements  within  the  element 
is  given  by 


U  =  <*,  *"  *1  *  +  V 
v  =  <*4+  *s  *  +  *<.  y 


(A-l) 


where  the  constants  are  determined  from  the  six  nodal  displacements  and 
nodal  coordinates  as 


and 


where 
and  c 


is  the  in-plane  area  of  the  element.  The  coefficients 
m 

are  given  by 


(A-2) 


(A-3) 


a  .  ,  b  . , 
1  l 


^4  “  Yk  - 


151 


(A-4) 


b;  =  Vj  -  yK 

C  i  “  *K  “ 

where  x  and  y  are  coordinates  of  the  nodal  points.  The  other  coefficients 
for  subscripts  "j"  and  "k"  are  obtained  by  cyclic  permutation  of  the 
subscripts  i,  j,  and  k. 

Element  Strain.  The  total  strains  at  any  point  within  an  element  are 
defined  in  terms  of  the  displacement  derivatives  as 


From  Eqs  A-l  to  A-5,  the  total  strains  are  written  in  terms  of  nodal  displacements 
and  coordinates  as 


it}--  t6M 


where  {U}  is  the  generalized  nodal  displacement  {Ur  =  (u.v.u.v.u,  v.) 

i  i  J  1  k  k' 


and 
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0 
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hi 

Ck  J 

The  superscript  T  denotes  the  matrix  transpose. 


ELASTIC  ANALYSIS 

For  linear  elastic  and  isotropic  materials,  the  relationship  between 

stresses  {a},  strains  (e } ,  initial  stresses  {a  1  and  anv  initial  strains 

o 

{e^}  is  given  by 


[°H  f<l- Wp  H] 


(A- 8) 
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where  [D  ]  is  the  elastic  material  property  matrix.  The  matrix  [D  ]  for 


plane-stress  conditions  where  =a  xz  =cr  ^£>0  is  given  bv 


1  V 

~Y  I 

0  O 


o  1 

o 


I- v 
a 


(A-9) 


where  E  and  v  are  the  modulus  of  elasticity  and  Poisson's  ratio,  respectively 
For  plane-strain  conditions  where  ez=  0,  the  elastic  material  property 
matrix  is  given  by 


L°] 


fl+v)(l 1-2  v) 


V 

//-l 


Vi-v 


O 

l-Zv 
Z(i-v)  J 


(A-l 0) 


Under  plane-strain  conditions,  a  normal  stress  also  exists  and  is  given  by 


r  1/  (<rx  -f  cry)  (A-H) 

Method  of  Solution.  The  equation  which  governs  the  elastic  response 
of  a  discretized  structure  can  be  derived  from  the  principle  of  virtual 
work  [31]  and  is  given  by 

Ik]  {u]  -  { p\ *  ld\ 

where  [K  ]  is  the  elastic  stiffness  matrix  of  the  structure,  I  is  the 
generalized  displacement  vector,  {Pi  is  the  externa]  applied  load  vector, 
and  {Q}  is  the  force  vector  due  to  the  presence  of  initial  stress  and  or 

initial  strain. 
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The  coefficients  of  the  elastic  stiffness  matrix  are  obtained  from 


M  =  I  /[b]T[o][b] 

frt  z  I 


(A-13) 


where  the  integration  is  taken  over  the  volume  of  each  element  and  the 
summation  is  over  all  elements  in  the  structure.  The  nodal  forces  due  to 
initial  stresses  are  given  by 

[q\  ~  2  (A_i4) 

hir  i 

and  the  nodal  forces  due  to  initial  strains  are  given  by 

° 

ELASTIC-PLASTIC  ANALYSIS 

The  application  of  the  finite  element  method  to  problems  involving 
materials  that  obey  linear  constitutive  laws  is  straight  forward  because 
the  material  properties  are  constant.  Therefore  only  one  solution  is 
required  to  obtain  displacements  for  the  elastic  structure.  However,  for 
elastic-plastic  problems  the  coefficients  in  the  stiffness  matrix  vary  as 
a  function  of  loading.  Thus,  the  elastic-plastic  displacements  are  usually 
obtained  by  applying  small  load  increments  to  the  structure  and  updating 
the  coefficients  of  the  stiffness  matrix.  Another  technique  called  the 
"residual  force"  method  [11]  avoids  modifications  of  the  stiffness  matrix 
by  adding  on  a  so-called  plastic  load  vector  to  the  force  side  of  the 
equilibrium  equation  (i.e.  Eq  A-12) .  Only  the  residual  force  method  will 
be  discussed  herein. 

Yield  Criterion.  In  any  elastic-plastic  material  the  elastic  formulation 
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can  be  used  prior  to  plastic  yielding.  Thereafter  it  is  necessary  to  have 
a  yield  criterion  to  determine  the  state  of  stress  at  which  yielding  occurs. 
The  von  Mises  yield  criterion  is  one  of  the  most  widely  used.  This  criterion 


assumes  yielding  is  caused  by  the  maximum  distortion  energy  [64].  The 
yield  criterion  for  plane  stress  conditions  is  given  bv 

F  =  F H  1 IV+  *-«  «t  +  3  -,;]t  ir 

and  for  plane  strain 

J. 

F=  "[<rxt+<ryl+0iZ"'rx<ry-ci^-  +  (a-17) 

where  a  is  the  uniaxial  yield  stress.  If  the  state  of  stress  is  such  that 
F  <  0,  the  material  is  still  in  the  elastic  range.  When  F  =  0,  a  plastic 
state  is  obtained  and  one  of  the  flow  theories  of  plasticity  must  be 
employed  to  determine  subsequent  plastic  behavior  under  increasing  stress 
or  strain. 

Flow  Theory.  One  of  the  basic  assumptions  in  the  theory  of  plasticity 
is  that  the  total  strain  {e }  or  total  strain  increment  {dc1'  can  be  decomposed 
into  elastic  and  plastic  strain  components  as  follows: 

or,  incrementally. 

In  the  incremental  theory  of  plasticity  the  plastic  strain  increment  vector 
{de*1}  is  a  function  of  the  current  state  of  stress  and  is  related  to  the 
yield  criterion  through  the  Theory  of  Plastic  Potential  [80] 

i^€?\  ~  ^  {  5H 

where  >,  is  a  positive  scalar  quantity.  This  flow  law  is  also  written  in 
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(A-18) 


(A- 19) 


\ 


V 


terms  of  strain  rate 


{cUP3r  l'  falTjl  (A-21) 

In  this  case  has  the  significance  of  the  coefficient  of  viscosity. 

Equations  (A-20)  and  (A-21)  are  also  known  as  Drucker's  Normality  Principle  [64] 
which  by  its  name  specifies  that  the  plastic-strain  increment  vector  is  to  be 
aligned  normal  to  the  yield  surface  in  nine  dimensional  stress  space.  When 
the  von  Mises  yield  criterion  is  used  with  Eq  A-20  the  resulting  expression 
for  (deP)  is  identical  to  that  proposed  by  Prandtl  and  Reuss  [64].  The 
total  strain  increment  vector  can  now  be  written  as 

where  the  elastic  strain  increment  vector  has  been  related  to  the  stress 
increments  {da}  through  the  elasticity  matrix.  Therefore,  if  2.  was  known, 
then  the  desired  stress-strain  relation  for  an  elastic-plastic  material 
would  be  obtained.  When  yielding  is  occurring,  the  total  differential  of 
Eq.  A-16  or  Eq.  A-17  gives 

I7 

dr-  -  c J  (r  =  o  (a-23) 

The  increment  in  yield  stress  da  is  obtained  from  a  uniaxial  tensile  test  as 

=  (j ^)dep=  H'Jer  <A-24) 

where  I!'  is  the  slope  of  the  stress-plastic  strain  curve,  daU  is  the 
uniaxial  stress  increment,  de^  is  the  uniaxial  plastic  strain  increment. 

Using  Eq .  A-20  for  the  uniaxial  case  gives  dr^  =2  ■  Thus.  I.q .  A-21  becomes 


157 


Eliminating  ^  from  Eqs  A-22  and  A-25  resn’*.,  j.n  an  explicit  expression 


relating  increments  in  stress  to  increments  in  total  strain  [44].  This 
expression  is 


-  lDEp] 


(A-26) 


where 


tDep]=  l°M°l 


(A- 2  7 ) 


The  matrix  [D  1  is  the  elastic-plastic  matrix  which  replaces  the  elasticity 
EP 

matrix  [D  ]  in  an  incremental  analysis.  For  an  elastic-perfectly  plastic 
material,  H '  is  set  equal  to  zero.  In  general  the  slope  of  the  uniaxial 
stress-plastic  strain  curve,  H',  varies  with  plastic  strain.  Therefore 
to  relate  a  mul tiaxiai  plastic  strain  state  to  a  uniaxial  experimental 
stress-plastic  strain  curve,  an  effective  plastic  strain  is  defined  in 
incremental  form  as 


ELASTIC-PLASTIC  SOLUTION  TECHNIQUES 

The  procedures,  used  to  solve  small  displacement  elastic-plastic 
problems  incrementally  within  a  finite  element  computer  program,  may  be 
divideci  into  two  categories.  In  one,  the  effects  of  plasticity  are 
accounted  for  directly  in  the  stiffness  matrix.  The  second  category  treats 
plastic  behavior  as  an  additional  plastic  load  that  is  combined  with 
applied  or  external  loads  in  the  equilibrium  equation  (i.e.  Eq .  A-12) . 

These  two  procedures  are  referred  to  as  the  "tangent  Modulus"  and  "residual 
force"  methods  respectively.  Only  the  residual  force  method  in  the  form 
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of  "initial  stress"  and  "initial  strain"  will  be  summarized  herein. 

Initial  Stress  Method  [44].  The  equation  which  governs  the  response 
of  a  discretized  structure  under  loads  which  cause  plastic  deformation  [3]]  is 


=  f  Pi '  +  £»lx' 


(A- 29) 


where  [« . ,]  is  the  elastic  stiffness  matrix,  {U1  is  the  generalized  displacement 
vector,  {P}  is  the  applied  load  vector,  and  {Q }  is  the  "effective"  plastic¬ 
load  vector  which  accounts  for  elements  in  a  plastic  state.  The  initial 
stress  method  approaches  the  solution  to  an  elastic  plastic  problem  by 
applying  a  series  of  small  load  increments  to  the  structure  until  the 
desired  load  is  reached  ({P}^  =  {P } ^  ^  +  {dP  t ) .  The  superscript  i  denotes 
the  current  increment  and  i-1  denotes  the  preceding  increment.  After  each 
load  increment  an  iterative  process  is  required  to  stabilize  the  plastic¬ 
load  vector.  The  subscript  I  denotes  the  current  iteration  and  1-1  denotes 
the  preceding  iteration.  During  the  ith  increment  a  purely  elastic  problem 
is  solved  and  the  increments  in  total  strain  { d e:  1  and  corresponding  elastic 
stress  {dc^}  are  computed  from  the  displacement  increments  {dll}  for  every 
element.  Because  of  the  material  nonlinearity  the  stress  increments  are 
not,  in  general  correct  or  if  the  correct  stress  increment  for  the  corresponding 
strain  increment  is  {do},  then  a  set  of  body  forces  or  plastic-load  vectors 
{ dQ }  caused  by  the  "initial"  stress  (dao )  =  ({da^l  -{dal)  is  required  to 
maintain  the  stress  components  on  the  yield  surface  or  compatible  with  the 
uniaxial  stress-strain  curve.  The  correct  stress  increment  is  computed 
with  Eq  A-26.  The  plastic  load  increments  are  computed  from  Eq  A-14 


(A- 30) 


Elements  are  in  the  elastic  state  when  {do  1  =  0.  The  total  plastic-load 
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vector  is  then  computed  as 


falj  =  + 


(A-31) 


At  the  second  stage  of  computation  the  new  force  system  {Q}^  is  added  to 
the  applied  load  vector  and  a  new  set  of  displacements  is  obtained. 

Again,  some  of  the  stresses  are  likely  to  exceed  the  yield  criterion  and  a 
new  set  of  plastic-load  increments  are  computed.  This  iteration  process 
is  repeated  until  the  change  in  the  plastic-load  vector  is  sufficiently  small. 
See  Fig  A-2  for  a  uniaxial  schematic  of  this  iterative  procedure  and 
Fig  A-3  for  the  mathematical  algorithm.  Consider  points  A,B,  and  C  in 
Fig  A-2.  Point  A  is  the  state  prior  to  the  load  increment.  Point  B  is 
the  state  after  the  load  increment  has  been  applied  and  one  "initial"  stress 
iteration  has  been  accomplished.  Point  C  is  the  state  of  stress  and  strain 
sought  after  which  satisfies  equilibrium  with  external  loads  and  compatabili tv 
with  the  material's  stress-strain  curve.  Notice  point  B  satisfies  compatabilit 
but  not  equilibrium  since  {dolj.  j<  {do^l^. 

INITIAL  STRAIN  METHOD 

The  initial  strain  method  parallels  the  initial  stress  method  somewhat 
and  accordingly  this  development  will  begin  just  after  step  4  in  the 
"Initial  Stress"  algorithm  in  Fig  A-3. 

The  elastic-plastic  material  matrix  [D  ]  is  used  as  follows 

EP 


(A-32) 


This  plastic  strain  increment  {dfP  :  is  then  used  to  calculate  a  plastic 
force  vector  increment 


Mr  I 


(A- 3  3) 
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Figure  A-2 . 


Uniaxial  "Initial"  Stress  Iteration  Schematic 


For  Each  Load  Increment  dP  in  Plastic  Range  44 
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Figure  A-3.  Initial  Stress  Algorithm 
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This  plastic  force  vector  increment  is  added  to  the  external  force  vector 
increment {dP }  for  the  augmented  global  force  vector  used  in  the  next 
iteration  as  follows: 


- 1 *Y't  fJ <*-«> 
{Jelx*,'  <A'35) 

The  stress  increment  {do}  is  calculated  as  follows 

<A-36) 

Steps  (1)  -  (4)  of  the  Initial  Stress  algorithm  and  Eqs  A-32  through  A-36 
above  are  repeated  until  compatability  with  the  materials  stress-strain 
curve  is  established.  Compatability  is  shown  to  be  achieved  after  "n" 
iterations  in  Fig  A-4.  Also,  compatability  would  display  itself  by 
little  or  no  change  in  the  plastic  strain  increment  between  iterations. 

Note  that  equilibrium  is  continually  satisfied  in  this  initial  strain  method. 
This  version  of  the  initial  strain  method  differs  from  Marcal  [45]  by  the 
fact  that  iterations  within  a  load  increment  are  done  globally  rather  than 
within  each  element  as  the  "constant  strain"  method  of  iteration  implies 
in  Marcal 's  paper. 


STRAIN 


A-4.  Uniaxial  "Initial"  Strain  Iteration  Schematic 


APPENDIX  B.  ITERATIVE  SOLUTION  TECHNIQUE 
FOR  NODE  POINT  DISPLACEMENTS 

The  matrix  equation  which  governs  the  response  of  a  discretized  structure 


is 


M  {u\  --  {?} 


(B-l) 


where  [  Kj  is  a  symmetric  positive-definite  n  x  n  matrix,  {U }  is  the  unknown 
node  point  displacement  vector,  and  {P}  is  a  known  load  vector.  In  the 
finite  element  method  for  structural  analysis,  the  matrix  [K  ]  is  usually 
highly  banded  and  if  stored  in  compacted  form  (i.e.  only  nonzero  terms 
retained)  requires  much  less  space  in  the  computer  than  the  product  n  x  n 
reflects.  Also,  if  there  are  changing  boundary  conditions,  such  as  freeing 
nodes  to  simulate  crack  growth,  then  the  [  |<J  matrix  must  be  recomputed. 

A  solution  technique  that  works  well  with  4compacted  [  kJ  matrices  and 
conveniently  admits  boundary  condition  changes,  is  the  Gauss-Seidel  iterative 
technique  with  over-relaxation  [14],  This  technique  may  be  implemented 
in  the  following  manner[fi7].  Consider  Eq  B-l  rewritten  as 


1*7  = 


5U  1  SXY' 

_ i _ _ 

SYX  !  SYY 


( B— 2 ) 


where  U^  and  U^  represent  node  point  displacement  vectors  in  the  x  and  y 
direction  respectively,  P^  and  P  represent  the  node  force  vectors  in  the 
x  and  y  direction  respectively.  The  submatrices  SXX,  SXY,  SYX,  and  SYY 
in  the  matrix  [  k]  have  dimensions  x  ^ ,  but  due  to  their  bandedness 
can  be  compacted  to  a  matrix  which  is  x  9.  The  dimension  9  minus  1 
reflects  how  many  adjacent  nodes  can  be  connected  to  any  given  node.  This 
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is  not  very  restrictive  since  triangular  finite  elements  develop  undesirable 
aspect  ratios  if  there  are  more  than  8  nodes  connected  to  any  one  node. 

Appropriate  terms  of  the  matrix  [  K]  are  retrieved  from  the  compacted 
submatrices  with  the  help  of  a  matrix  NP  and  a  vector  NAP.  The  vector  NAP 
has  the  dimension  ^  and  the  matrix  NP  has  the  dimension  ~  x  9 .  The  Ith 
component  of  NAP  or  NAP(I)  stores  the  number  of  adjacent  node  points 
connected  to  nodepoint  I.  The  (I,  J)  component  of  NP  or  NP(I,J)  stores  the 
address  of  the  terms  in  the  submatrices  associated  with  the  Jth  adjacent 
node  point  connected  to  node  point  I.  Note  that  for  node  I,  J  may  go  from 
1  to  NAP(I) . 


Consider  the  governing  equation  for  node  point  1  displacements  written  as 


\  sxx  a,  0  sxr(i,0 

SYXCX,  0  SYY(V) 


VA  P  (l) 


[sxx(x,j)  sxr(x,i?  i 

SYX(l,3)  SVY  Ci, 3) 


If  the  right  hand  side  of  Eq  B-3  is  defined  as  the  ver  tor }  t  hen  solving 
for  the  displacements  at  node  I  yields 


UyC*) 


(  >s-  (  I 


1 


dtCx)J 

Uy£r)f 


SXX  (1,0  SXY(I,|) 
SYx  (X, »)  SYY(I,0. 
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(B-4) 


Note  that  the  matrix  to  be  inverted  in  Eq  B-4  is  only  a  2  x  2.  Also  since 
this  is  the  only  place  these  terms  of  the  submatrices  are  used  this  2x2 
may  be  inverted  and  its  components  stored  in  their  original  submatrix  locations. 
To  incorporate  an  over-relaxation  factor,  Eq  B-4  is  modified 
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as  follows 


tsxxdv)  sxY(r,i) 

SYX(X,|)  SYY(I;|) 


(B-5) 


where  the  superscripts  m  and  m-1  refer  to  iteration  number.  The  left 
hand  side  of  Eq  B-5  is  the  change  in  displacements  between  iterations 
without  applying  an  over-relaxation  factor.  But  the  new  total  displacements 
for  iteration  m  using  an  over-relaxation  factor  are 


m-l 

+  XFAC 


A  U%  (l) 
A  Ui  Cx) 


(B-6) 


where  XFAC  is  the  over-relaxation  factor  which  normally  ranges  from  1.8  to 
1.9  for  structural  analysis. 

Convergence  of  these  iterations  is  checked  bv  computing  an  effective 
force  unbalance  term,  SUM,  defined  as 

nJz 

SUM  =  2.  [(  AUx£i)  *SXX  (i}  i)(  +  (Atfy60*srrcr,|)|J  (B'7) 

i=( 

If  SUM  becomes  less  than  a  specified  small  value,  9,  iterations  are  stopped 
and  the  node  point  displacement  solution  is  obtained.  The  value  of  9  is 
chosen  based  on  examining  solutions  for  various  sizes  of  9.  A  good  starting 
value  for  9  is  one  tenth  of  the  applied  load.  The  final  9  is  then  chosen 
based  on  the  amount  of  accuracy  desired. 

Displacement  boundary  conditions  are  easily  input  to  this  solution 
routine  by  simply  specifying  the  desired  node  point  displacements  and  then 
require  the  iteration  algorithm  to  skip  ove-  the  node  displacement  equations 
which  have  fixed  displacements.  Likewise,  if  a  fixed  node  is  released 
during  the  analysis  such  as  for  modelling  crack  growth,  this  node's  dis- 
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placement  equation  may  be  reactivated  in  the  iteiation  algorithm. 

Convergence  of  this  solution  technique  (i.e.  when  SUM  <_ 6)  is 
dependent  on  the  initial  guess  for  the  node  displacements.  Usually  for 
convenience  all  unknown  displacements  are  initialized  at  zero.  However, 
for  each  succeeding  solution,  such  as  in  a  nonlinear  incremental  analysis, 
much  better  initial  displacement  values  are  available  from  the  prior  solution. 
These  initial  displacements  from  prior  solutions  significantly  reduce  the 
number  of  iterations  to  convergence. 
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APPENDIX  C.  THE  J  AND  C*  INTEGRALS 


J  INTEGRAL 


Rice's  J  Integral  is  defined  as 


where 


T=  J"[vV (e£Jr  -7]  Js] 


A  ^rrrn 

W(f.j')-  4  Oijde.j 


(C-I ) 


(C-2) 


and  r  is  a  closed  loop  around  the  cracktip  as  shown  in  Fig  II-2.  Expanding 
J  and  integrating  along  a  rectilinear  path  in  the  x  and  v  directions  results  in 


j  «  Jt  <W  -  n  §»  -  <r  *)Jy  -  ft +  ^SL)  Jx  1 


(  C—  3) 


The  following  describes  a  numerical  procedure  for  calculating  the  J  integral 
with  a  finite  element  program  that  incorporates  constant  strain  triangles. 
Consider  Fig  C-l  which  is  a  region  of  elements  taken  from  a  finite  element 
model  of  a  cracked  geometry.  Paths  1  and  2  are  two  possible  paths.  The 
contribution  to  path  1  from  element  2  is 


a  i 1  -  <5  « I?  <*>  -  f  r  <*>  It  f*)3  -  y 


(C-4) 


where  the  number  inside  parenthesis  refers  to  the  element  number  these 
values  came  from.  A  similar  contribution  would  come  from  element  5  for 
path  1.  Path  2  which  also  runs  through  element  2  would  also  have  the  above 
contribution  from  element  2  but  as  path  2  turns  and  runs  along  the  border  of 
elements  6  and  1  the  following  is  the  contribution 


+  ^y0)^0)^-y0)^0)]  (~^) 
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Notice  in  this  case  an  average  of  the  stresses  and  strains  in  the  two 


elements  is  taken  by  effectively  running  half  the  element  length  in  element 
6  and  half  in  element  1. 

The  strain  energy  term  U(2)  in  Eq  C-4  is  calculated  as  follows 

c-  o)  e*  o)  +  J*  cT"  (c-6) 

S  4  o  d 

Although  these  equations  are  shown  for  specific  element  numbers  their  form 
is  used  for  all  elements  along  the  J  integral  path. 

C*  INTEGRAL 

The  C*  Integral  [23]  defined  as 


C*  ^  -  T;  iili] 


(C-7) 


where 


;p 


w 


/■s  ^  nn 

-  S  J  e,F 

0  'J  ' 


(C-S) 


and  P  is  the  same  type  of  path  as  for  the  J  integral.  C*  may  be  obtained 
by  replacing  strain  and  displacement  in  the  J  integral  with  strain  rate 
and  displacement  rate  respectively.  The  rational  behind  this  is  based  on 
the  assumption  that  the  material  being  analyzed  behaves  as  a  creeping 
solid  such  that 


*  f* 

€  = 


X  (<r) 


(C-9) 


or  in  a  multiaxial  sense 


3  w 


(C.O 


(C-10) 
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where  W*  was  defined  in  eq .  C-8.  Notice  Eq.  C-9 ,  which  is  considered  a 


plastic  strain  rate,  makes  no  provision  for  admitting  elastic  strain  rates. 
Hence  the  claim  of  path  independence  for  the  C*  integral  can  only  be 
approached  in  a  realistic  elastic-plastic  material  once  the  stresses  have 
reached  a  steady  state  value  (i.e.  =  0).  With  this  restriction  on 

stresses,  C*  may  also  be  defined  as  J  or  the  time  rate  of  change  of  the 
J  integral.  In  general  J  would  include  the  time  derivatives  of  stress  and 
traction  but  if  they  are  restricted  to  zero  then  j  is  equal  to  C*  for  a 
creeping  solid.  Also  with  a  equal  to  zero,  Eq  C-8  may  be  directly 
integrated  to 


W*  =  (C-11) 

Using  this  form  of  W*,  C*  may  be  numerically  calculated  in  the  same  fashion 
as  the  J  integral  with  the  exception  of  using  strain  rates  and  displacement. 
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APPENDIX  D.  DETERMINATION  OF  BODNER 


MATERIAL  MODEL  CONSTANTS 

This  appendix  describes  how  each  of  the  material  constants  for  the 
Bodner-Partom  material  model  can  be  determined.  Normally  high  costs 
and  material  shortages  prevent  obtaining  more  than  one  or  two  stress- 
strain  curves  and  the  same  number  of  creep  tests.  Ideally,  to  best 
characterize  the  material,  several  stress-strain  curves  should  be 
generated  over  a  wide  range  of  strain  rates,  similar  to  Fig  III-2. 
Likewise,  several  creep  tests  should  be  performed  over  a  wide  spectrum 
of  stress  levels. 

In  general  the  Bodner  material  constants  are  dependent  on  temperature. 
However,  the  temperature  dependence  is  suppressed  by  performing  the 
material  characterization  tests  (i.e.  stress-strain  and  creep)  at  the 
same  temperature  that  the  Bodner  model  will  be  applied. 

To  determine  the  Bodner  constants  from  uniaxial  test  data,  the 
Bodner  equations  are  written  in  uniaxial  form  as  follow:  (the  total 
strain  rate  is  the  sum  of  plastic  and  elastic  strain  rates) 

6  -  «f  €?  (D-i) 

and  the  plastic  strain  rate  for  uniaxial  tension  is 


(D-2) 


where 


?r  1,  -  -  ?c)  «xp[-m  VN/p]  (D-3) 


(D-5) 


ir..  -  -A  (¥)'., 

-3 

over  small  time  intervals,  such  as  a  typical  stress-strain  test  at  10  sec 
strain  rate,  the  recovery  term,  Zrec,  may  be  ignored  and  then 

*  •  p 

Wp  =  O'  £  (D-6) 

Since  test  data  can  be  resolved  into  the  forms,  and  a,  Eq.  E-2 

is  solved  for  Z  which  then  is  a  function  of  t  and  o  as  follows 

P 

<r[;rrr  '"  (%■&■)]*"  (Mi 

The  viscoplastic  material  constants  in  these  equations  are  broken 
into  "short  time  response"  and  "creep"  groupings  for  determination. 

SHORT  TIME  RESPONSE  CONSTANTS 

The  short  time  response  constants  for  the  Bodner  model  are  Dq,  n, 
m,  Zq,  Z^.  These  constants  are  primarily  determined  by  using  stress- 
strain  test  data. 

4  -1 

The  constant  Dq  is  normally  assigned  the  value  of  10  sec  .  For 
high  strain  rate  applications,  Dq  may  be  set  higher  (e.g.  106sec-1) 
which  would  result  in  small  changes  to  the  other  constants. 

The  constant  n  is  directly  related  to  the  model's  strain  rate 
sensitivity.  High  n  values  reflect  low  strain  rate  sensitivity  and 
vice  versa.  Changes  in  n  affect  the  stress  values  for  a  given  strain 
rate  by  shifting  up  or  down  the  family  of  stress-strain  curves,  but  the 
shape  of  the  curve  is  preserved.  The  value  of  n  is  determined  in  an 
Iterative  fashion.  The  first  estimate  of  n  should  be  between  1  and  10 
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based  on  past  material  modeling  efforts  [8].  Plots  of  Z  versus 

are  then  made  by  simultaneous  use  of  each  stress-strain  curve,  Eq  E-l, 

Eq  E-6,  and  Eq  E-7.  The  value  of  n  is  then  adjusted  with  the  objective 
of  making  all  Z  versus  curves  from  each  stress-strain  test  fall  on 
top  one  another.  This  value  of  n  will  then  satisfy  the  requirement 
that  Z  is  a  single  value  function  of  as  given  in  Eq  E-3. 

The  first  approximation  of  Zq  can  also  be  determined  from  Eq  E-7. 

A  small  value  of  (e.g.  10  ^sec  *)  and  the  lowest  apparent  initial 
yield  stress  from  the  stress-strain  test  data  are  substituted  into  Eq  E-7. 
The  resulting  value  of  Z  is  defined  as  Zq.  Note  that  Zq  is  the  primary 
constant  that  determines  the  stress  level  at  which  significant  plastic 
straining  (i.e.  d*5  10  ^sec  *)  begins. 

The  constants  Z^  and  m  are  determined  by  rewriting  Eq.  E-3  as 

In  (*x~  =  In  (ix-  20)  -  ^  Wp  (D-8) 

An  iterative  process  is  now  begun  to  determine  Z^.  The  first  estimate 
of  Z,  should  be  larger  than  Z  (e.g.  1.5  Z  )  since  Z,  is  the  maximum 
value  for  Z.  By  incorporating  this  estimate  of  Z^,  a  plot  of  In  (Zj-Z) 
versus  Wp  is  made  based  again  on  stress-strain  test  data,  Eq  E-l,  E-6, 
and  E-7.  This  should  approximate  a  straight  line  whose  slope 
is  the  constant  m  and  extrapolation  of  the  line  to  ■  0  provides  a 
value  of  Zq.  If  this  ZQ  obtained  graphically  does  not  agree  with  the 
previous  value  for  Zq,  adjust  Z^  accordingly  and  reiterate. 

These  values  for  n,  m,  Zq,  and  Z^  which  primarily  govern  the 
short  term  stress-strain  behavior  should  be  input  to  a  computer  program 
that  numerically  integrates  Eqs  E-l  through  E-4.  Stress-strain 


predictions  from  the  model  should  be  made  for  each  experimental  strain 
rate  to  see  how  good  the  fit  is.  Only  minor  adjustments  to  the  short-time- 
response  constants  should  be  necessary  to  best  fit  test  data. 

CREEP  CONSTANTS 

The  creep  constants  for  the  Bodner  model  are  Z i.  A,  and  r.  These 
constants  are  determined  based  on  data  from  at  least  two  creep  tests  at 
different  stress  levels. 

During  secondary  creep,  when  plastic  strain  rate  is  approximately 

constant,  it  follows  from  Eq  E-7  that  Z  must  also  be  a  constant.  In 

addition,  if  Z  is  a  constant  it  follows  from  Eq  E-3  that  W  is  a  constant 

P 

which  makes  *  0.  Hence, with  =  0  in  Eq  E-4  and  combining  with 

Eq  E-5 

<r«r  =  A(Hr)  <D",) 

The  constant  Z^  represents  the  minimum  value  for  Z  corresponding 
to  secondary  creep.  If  creep  occurs  below  the  apparent  yield  stress 
implied  by  ZQ,  the  value  of  Z^  must  be  less  than  Zq.  Moreover,  Z^ 
must  be  less  than  or  equal  to  the  smallest  value  of  Z  from  Eq  E-7 
when  using  creep  test  data  (i.e.  and  0).  After  selecting  a  value 
for  Z^,  the  constants  A  and  r  can  be  determined  after  rewriting  Eq  E-9 
in  terms  of  natural  logarithms 

In  (r€r)+  *  In  (A  +  r  (d-io> 

Stress  and  plastic  strain  rate  are  substituted  into  Eq  E-10  along 


m 


with  the  appropriate  Z from  Eq  E-7  and  creep  test  data.  With  data  from 
creep  tests  at  two  different  stress  levels  two  linear  equations  in  terms 
of  the  two  unknowns  A  and  r  are  developed.  The  constants  A  and  r  are 
then  determined  by  simultaneous  solution  of  these  two  equations  based  on 
Eq  E-10  and  creep  test  data. 

The  complete  Bodner  model  should  then  be  tried  out  in  a  computerized 
numerical  Integration  scheme  in  an  effort  to  make  final  adjustments 
to  the  constants  for  a  best  fit  to  the  test  data. 

It  became  apparent  to  the  author  that  in  working  with  this  material  model 
a  sensitivity  study  is  required  in  order  to  see  the  effect  of  changes  in  the 
respective  constants  on  creep  crack  growth.  Yet  it  is  felt  that  the  results 
and  conclusions  presented  in  the  main  body  of  this  report  provide  realistic 
directions  and  trends  of  creep  crack  growth. 
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